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CHAPTER  1 


INTRODUCTION 


1.1  Introduction 

Texture  is  important  characteristic  for  the  analysis 
of  many  types  of  images.  It  is  an  important  feature  for 
discrimination  and  identification  of  regions  in  images  and 
as  a  result  the  vast  majority  of  work  on  texture  has 
concentrated  on  these  applications.  Although  many 
different  texture  discrimination  techniques  have  been 
developed,  most  are  ad  hoc . 

The  problem  inverse  to  texture  analysis  is  texture 
synthesis,  or  the  generation  of  image  fields  having 
analytical  and  visual  characteristics  similar  to  natural 
textures.  Texture  synthesis  has  been  over-shadowed  by  the 
emphasis  placed  on  the  discrimination  problem  and  its 
applications.  Little  work  has  been  done  on  synthesis  even 
though  numerous  applications  exist.  For  example, 
intelligent  image  sensors  could  transmit  boundaries  of 
textured  image  regions.  Based  on  statistics  gathered  by 
the  sensor,  this  region  could  be  reconstructed  using 
simulation  techniques  with  little  or  no  loss  of 


1 


information.  The  result  is  excellent  data  compression. 


Texture  synthesis  can  also  be  used  as  a  texture 
analysis  tool  leading  to  a  better  understanding  of 
textures  and  their  perception  by  humans  as  v/ell  as 
improved  methods  of  discrimination.  By  carefully 
controlling  the  statistics  of  a  texture  in  a  synthesis 
process  visual  changes  in  texture  are  observed.  Thus, 
texture  synthesis  methods  allow  researchers  to  identify 
and  measure  the  information  content  of  individual 
statistical  measurements.  By  assembling  these 
measurements  and  incorporating  them  into  a  texture 
simulation  process,  statistics  may  be  measured  from  a 
parent  texture  and  used  to  produce  a  texture  simulation. 
The  degree  to  which  the  parent  and  simulation  are  visually 
similar  indicates  the  value  of  the  statistical 
measurements  and  the  model  used  in  the  simulation  process. 
Given  a  group  of  statistical  measurements  which  are 
proposed  to  be  useful  texture  measures,  possibly  the  best 
may  be  chosen  based  on  the  quality  of  the  corresponding 
texture  simulations.  In  this  way,  researchers  are  able  to 
develop  better  d iscr im inat ion  as  well  as  better  synthesis 
methods. 
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1.2  Concepts  of  Texture  Synthesis 

Despite  its  importance,  a  precise  definition  of 
texture  does  not  exist.  Texture  is  often  considered  to  be 
composed  of  a  set  of  primitives  and  their  spatial 
organization.  More  important,  texture  usually  possesses 
an  invariance  property.  We  will  use  this  invariance 
property  as  one  definition  of  texture  in  this  thesis.  An 
observer  should  detect  no  visual  difference  between  one 
windowed  portion  of  a  textured  region  and  another.  Thus 
texture  is  also  a  function  of  window  size.  If  a 
difference  over  a  region  is  detected  then  either  the 
texture  is  not  homogeneous  (see  section  2.1)  or  a  larger 
or  smaller  window  should  be  used.  Windowing  is  very 
important  when  gathering  statistics  to  be  used  for  texture 
discrimination  or  texture  synthesis. 

The  approach  to  texture  synthesis  used  in  this  thesis 
is  outlined  in  Fig.  1.1.  As  a  first  step  in  the  synthesis 
process,  statistics  are  calculated  from  measurements  taken 
on  a  parent  sample  texture.  The  statistics  are  then  used 
to  estimate  model  parameters.  In  the  final  step,  these 
model  parameter  estimates  are  used  to  generate  a  texture 
synthesis . 

All  of  the  digital  images  in  this  thesis  are  512  by 


3 


512  pixels.  They  have  either  256  gray  levels  (continuous 
tone)  or  2  gray  levels  (binary) .  The  original  parent 
texture  images  in  this  thesis  have  been  cnosen  from  an 
album  by  Brodatz  [1],  High  quality  prints  obtained  from 
the  photographer  were  scanned  and  digitized  at  the  'JSC 
Image  Processing  Institute. 

The  performance  measure  for  any  texture  synthesis 
method  is  purely  visual.  Thus  the  rating  of  any  synthesis 
system  must  be  made  by  an  observer  and  is  subject  to  human 
variability.  These  facts  must  be  remembered  when 
considering  the  results  of  synthesis  techniques  discussed 
in  this  thesis  and  in  other  works.  For  this  reason,  the 
visual  analysis  of  results  in  this  thesis  is  left 
primarily  to  the  reader.  Nevertheless,  general  guidelines 
and  trade-offs  involved  in  texture  synthesis  have  been 
developed  in  the  course  of  this  work. 

The  success  of  any  texture  synthesis  as  well  as  any 
image  processing  technique  (enhancement,  restoration, 
etc.)  also  depends  on  the  display  medium  used  for  final 
results.  It  is  unfortunate  that  the  product  of  so  much 
work  is  subjected  to  the  imperfections  of  recording  and 
printing  processes.  We  have  attempted  to  minimize  these 
degradations  as  much  as  possible. 
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1.3  The  Stochastic  Model 


The  approach  used  to  synthesize  textures  in  this 
thesis  is  probabilistic  (statistical)  rather  than 
structural  because  the  structural  approach  would  probably 
require  a  different  processing  method  for  each  texture. 
Textures  which  are  irregular  and  highly  random  are  often 
difficult  to  analyze  using  a  structural  approach  as  they 
do  not  have  "structure"  as  such.  On  the  other  hand,  such 
textures  are  easily  analyzed  using  a  probabilistic 
approach.  Textures  which  are  not  composed  of  well-defined 
non-varying  primitives  usually  do  not  lend  themselves  to  a 
structural  approach.  Highly  regular  textures  which  are 
usually  best  explained  using  a  structural  approach  may 
also  be  studied  using  a  statistical  approach.  In  this 
case,  the  statistical  model  must  be  constructed  to  explain 
regular  and  periodic  events  with  less  importance  placed  on 
random  elements  of  the  model.  Naturally,  the  structural 
and  probabilistic  approaches  overlap  considerably. 

In  our  statistical  approach  to  texture  analysis  and 
synthesis  we  will  use  various  stochastic  models.  Textures 
are  analyzed  as  series  of  pixels  to  which  a  model  is  fit. 
The  time  or  space  series  of  successive  pixels  which  forms 
our  sample  parent  texture  is  regarded  as  a  sample 
realization  from  an  infinite  population  of  such  textures. 
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A  model  is  derived  which  attempts  to  explain  the  numerical 
sequence  of  pixels  in  the  observed  parent  texture.  This 
stochastic  model  may  then  be  used  to  simulate  a  similar 
sequence  which  becomes  our  texture  synthesis. 

The  power  of  the  stochastic  models  presented  in  this 
thesis  is  that  it  is  easy  to  use  the  model  in  a  mode  which 
synthesizes  textures  given  the  necessary  model  parameters. 
In  this  sense,  the  stochastic  approach  is  sufficient  to 
capture  everything  about  a  texture.  The  quality  of  the 
synthesis  is  also  a  measure  of  the  amount  of  information 
contained  in  the  model. 

1.4  Organization  and  Contributions  of  the  Thesis 

In  the  next  chapter,  a  one-dimensional  texture 
synthesis  model  is  mathemat ical ly  developed.  It  is 
analyzed  as  a  Markov  model  where  the  state  of  the  system 
is  considered  to  be  defined  by  the  sequence  of  previous 
pixel  values.  Similar  texture  synthesis  has  also  been 
done  by  Julesz  [8],  Purks  and  Richards  [3],  Conner s  [7], 
Abend  [6],  and  Gagalowicz  [4,5],  With  the  simple  model 
presented  in  Chapter  2,  textures  with  controlled 
statistical  properties  are  generated  and  used  to  study 
human  perception  of  texture. 

In  Chapter  3,  the  one-dimensional  model  is  extended 


7 


to  tw.  imensions  and  is  used  to  synthesize  natural  binary 
textures.  This  work  was  presented  earlier  by  Garber  [9] 
and  similar  work  was  done  by  Lu  and  Fu  [10],  This 
two-dimensional  probabilistic  model  which  generates 
textures  based  on  higher-order  probability  densities  is 
then  reduced  carefully  to  a  linear  autoregressive  model. 
Results  of  both  methods  applied  to  a  wide  variety  of 
textures  are  shown. 

In  Chapter  4,  a  method  to  reconstruct  higher-order 
densities  for  texture  synthesis  from  second-order 
measurements  is  presented.  The  results  confirm  the  value 
of  second-order  statistics  in  texture  analysis. 

In  Chapter  5,  the  autoregressive  model  is  applied  to 
continuous- tone  texture  synthesis.  Much  simpler  synthesis 
models  were  used  by  McCormick  and  Jaya ramamur thy  [2]  and 
Tou,  Kao  and  Cheng  [11].  The  success  of  their  method  was 
not  established  due  to  the  limited  number  of  textures 
involved  in  their  studies,  however  their  results  suggested 
that  time  series  models  could  be  used  to  generate  some 
natural  textures.  The  large,  two-dimensional,  linear 
model  presented  in  Chapter  5  is  then  extended  to  a 
quadratic  autoregressive  form  and  synthesis  results  with 
both  stationary  and  non-sta t iona ry  noise  are  presented. 
The  generated  textures  show  that  these  models  define 
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valuable  synthesis  methods. 

In  Chapter  5,  texture  synthesis  methods  wnich 
incorporate  more  than  one  autoregressive  model  are 
presented.  The  multiple-model  methods  methods  could  be 
valuable  in  synthesizing  textures  which  have  very  coarse 
structure  and  textures  which  are  composed  of  subtextures. 

In  Chapter  7 ,  a  texture  synthesis  method  which 
simulates  complete  probability  density  information  for  use 


in  the  synthesis 

process  is 

developed . 

This  method 

is 

computationally 

burdensome 

but  also 

yields  the 

best 

synthesis  results 

.  A  method 

similar 

to  this  will 

be 

valuable  as  processor  speed  increases  in  the  future. 

In  Chapter  8,  methods  for  adding  and  removing 
non-homogeneities  in  texture  mean  and  variance  are 
presented  . 

In  Chapter  9,  the  measures  used  to  estimate 
parameters  in  the  autoregressive  model  are  applied  to  the 
problem  of  texture  discrimination.  The  methods  differ 
from  those  presented  earlier  by  Deguchi  and 
Morishita  [12],  Kaiser  [13],  Pratt,  Faugeras  and 
Gagalowicz  [14],  and  Haralick  ejt  a  1 .  [15].  The  chapter 
illustrates  two  approaches  to  use  statistics  from  a 
synthesis  model  to  discriminate  textures. 
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The  experimental  results  of  this  thesis  are  largely 
visual  therefore  particular  attention  should  be  devoted  to 
the  figures.  A  casual  reader  could  read  section  2.1, 
Chapter  3,  Chapter  5,  and  Chapter  7  and  still  understand 
much  of  the  basic  concepts  of  the  work  as  well  as  major 
results.  Chapter  10  contains  a  final  summary  and 
comparison  of  the  models  and  their  corresponding  results. 

As  a  whole,  this  thesis  presents  results  in  texture 
simulation  using  methods  developed  herein  or  only  briefly 
mentioned  in  other  previous  studies.  The  texture 
syntheses  are  exceptional  in  some  cases  and  certainly 
notable  in  others.  This  work  should  encourage  additional 
research  in  the  field  of  texture  synthesis. 

1.5  Notation 

Any  variables  which  have  different  meanings  within 
this  thesis  are  clearly  defined  in  the  places  where  they 
are  used.  The  term  "complex"  means  complicated  rather 
than  a  variable  with  real  and  imaginary  parts.  The  terms 
"normal  distribution"  and  "normally  distributed"  are  to  be 
taken  in  a  statistical  probability  density  function 
(ie. Gaussian  distribution)  sense. 
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CHAPTER  2 


ONE-DIMENSIONAL  BINARY  TEXTURE  MODEL 

2.1  Introduction 

Texture  is  a  complex  image  attribute  that  has  been 
the  subject  of  much  research  and  is  difficult  to  define 
precisely.  The  relationship  between  discrimination  of 
textures  by  human  observers  and  the  mathematical 
attributes  of  textures  has  also  been  extensively 
researched.  Models  for  computer  discrimination  have  been 
proposed  based  on  statistical  parameters  considered  in 
some  aspects  to  be  primary  texture  measures. 

The  terminologies  in  a  portion  of  previous  texture 
work  have  often  been  vague  at  best.  As  a  result,  the 
terms  second-order  and  third-order  have  been  seriously 
twisted  and  misinterpreted  from  study  to  study.  In  this 
and  following  chapters,  we  will  attempt  to  suppress  this 
confusion  by  carefully  defining  the  various  terms. 

The  stochastic  approach  toward  texture  analysis 
considers  texture  fields  as  samples  of  two-dimensional 
stochastic  fields.  Assuming  that  we  are  dealing  with 
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sampled  and  quantized  imagery,  let  I(n^,n^2)  denote  the 
random  field.  Here  n^  and  n.^  are  integers  representing 
coordinates  of  points  in  the  image  plane.  Let  be  the 
vector  having  coordinates  n^  and  n^2  (  i  .e  .r^=  (n.^  ,n^2)  )  . 
Second-order  statistics  are  given  by  the  set  of 
second-order  joint  density  functions 


(Vf'Vj) 


(2.1) 


for  all  possible  vectors  n^  and  n j ,  where  and  Vj  are 
the  values  of  the  random  variables  I(n^)  and  I(rij), 
respectively.  In  most  texture  work  and  in  all  of  the  work 
in  this  thesis  (except  for  the  work  in  Chapter  7  and 
Chapter  8)  the  random  field  is  assumed  to  be  homogeneous, 
that  is,  all  orders  of  probability  densities  are  invariant 
through  translations.  Thus, 


P->-  -*■  =  P+  ( ?  ?) 

n.  ,n  .  n.+c,n.+c 

l  ]  I  D 

where  c  is  an  arbitrary  vector  constant.  As  an  example, 


p(v1,v2)  =  P(v3,v4) 


(2.3) 


where  V^,  V2,  V3  and  V4  are  as  shown  in  Figure  2.1.  In 
most  of  our  work,  dummy  values  of  random  variables 
(denoted  for  example  by  V^)  will  be  used  to  label  pixels 
at  vector  location  n;  . 
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Given  the  assumptions  that  a  texture  field  is 

homogeneous,  the  joint  density  functions  P  for  all  vector 

r 

separations  r  =  n^-nj  represent  the  most  complete  set  of 
second-order  statistics  possible.  The  statistical 
expectation  of  any  functions  of  these  joint  density 
functions  are  called  second-order  statistics.  If  a  pixel 
is  connected  to  any  of  its  neighbors  on  the  same  row,  that 
is,  if  we  consider  neighbors  immediately  to  the  left  or 
right  (such  as  Vg  and  Vg  in  Figure  2.1),  then  their  joint 
density  is  called  a  second-order  nearest-neighbor  joint 
density  and  any  statistical  expectations  of  the  joint 
density  are  second-order  nearest-neighbor  statistics. 
Nearest-neighbor  densities  and  nearest-neighbor  statistics 
are  very  important  in  this  chapter  as  the  textures  to  be 
generated  are  primarily  one-dimensional. 

Similarly,  third-order  statistics  are  given  by  the 
set  of  third-order  density  functions 


P-*-  -*■  -*•  (V.,V.,V.) 

ni'nj'nk  1  3  k 


(2.4) 


Assuming  homogeneity  of  tha  texture,  then 


n .  , n  .  ,n  n.+c,n.+c,n  +c  U.b) 

X  J  K  X  J  iv 

for  all  rti  and  an  arbitrary  vector  constant  <?.  As  an 
example , 
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in  Figure  2.2.  The  statistical  expectations  of  any 
function  of  these  third-order  densities  are  called 
third-order  statistics.  All  second-order  statistics  may 
be  derived  from  third-order  joint  densities.  In  this 
chapter,  third-order  statistics  involving  adjacent  pixels 
along  an  image  row,  such  as  the  pixels,  V^,  Vg,  Vg,  in 
Figure  2.2,  will  be  called  third-order  nearest-neighbor 
statistics . 

Julesz  [8]  created  computer  generated  patterns  with 
controlled  high-order  statistical  properties.  A 
conclusion,  often  referred  to  as  the  Julesz  conjecture, 
drawn  from  his  work  is  that  texture  fields  differing  only 
in  third-  and  higher-order  statistics  cannot  be 
discriminated  by  a  human  viewer.  Pollack  [16]  showed 
later  that  textures  whose  first-  and  second-order 
nearest-neighbor  probabilities  are  equal  may  be 
discriminated  by  varying  the  third-order  nearest-neighbor 
probabilities.  Purks  and  Richards  [3]  extended  this 
concept  to  create  texture  patterns  that  differ  only  in 
their  statistics  for  four  adjacent  points.  This  study 
indicates  that  such  textures  can  also  be  easily 
discriminated.  However,  as  was  pointed  out  by  Pratt  [14], 
the  second-order  probability  densities  of  the  two  fields 


are  not  constrained  to  be  equal  for  arbitrary  pixel  pairs 
along  an  image  line.  Thus  there  is  still  some  question  as 
to  the  relationships  between  measured  mathematical 
parameters  and  human  d  iscr iminabil ity.  Later  work  by 
Gagalowicz  [5]  seems  to  indicate  that  carefully  generated 
binary  patterns  whose  second-order  probability  pairs  are 
equal  for  arbitrary  distances  can  be  visually 
discriminated  by  human  observers  and  therefore  presents  a 
valid  contradiction  to  the  Julesz  conjecture.  Controlling 
different  statistics  of  a  texture  is  often  a  painfully 
difficult  process  and  as  a  result,  most  of  these  textures 
are  generated  using  approaches  unlike  the  Markov  approach 
of  this  chapter.  Often  blocks  of  pixels  are  generated 
with  certain  properties  or  patterns  with  special 
orientation  and  separation  are  used  to  study  the  effect  of 
statistical  changes  on  the  human  visual  system.  The 
mathematical  relationships  between  the  joint  density 
functions  of  any  texture  are  indeed  complex. 

For  these  reasons,  we  begin  with  one  rather  simple 
method  of  generating  Markov  one-dimensional  binary 
textures.  In  later  chapters,  these  ideas  will  be  modified 
and  extended  to  generate  and  simulate  two-dimensional 
textures. 

We  have  studied  in  detail  the  mathematical 


15 
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involved 


relationships  of  parameters  involved  in  binary 
computer-generated  one-dimensional  texture  patterns. 
Texture  patterns  in  this  paper  have  been  generated  using 
the  mathematical  relationships  derived  herein.  Methods 
have  been  developed  to  control  texture  statistics  for  both 
nearest-neighbor  and  non-nearest-neighbor  cases.  Examples 
of  both  types  of  textures  are  presented.  Using  these 
methods,  numerous  counter  examples  to  Julesz's  conjecture 
may  be  generated  and  are  illustrated  in  this  Chapter. 


Throughout  this  chapter,  the  P  (vi , V2 » • • • , VN )  will 
denote  the  nearest-neighbor.  Nth-order  joint  densities  of 
our  one-dimensional  texture  and  the  pixels  vi»v2'***,vN 
will  be  adjacent  to  one  another  in  the  sequence  as  shown 
in  Fig.  2.3  unless  otherwise  stated. 

2.2  Generation  Procedure 

One-dimensional  binary  textures  represent  the 
simplest  form  of  texture  possible.  It  is  believed  that 
such  binary  patterns  force  human  observers  to  utilize 
primitive  visual  mechanisms  in  discrimination.  They  are 
not  designed  to  replace  or  imitate  natural  textures  but 
are  experimentally  valuable  in  deriving  concepts 


concerning  texture  attributes  due  to  their  mathematical 
simpl icity. 


Figure  2.1  Second-order 
Statistics 


Figure  2.2  Third-order 
Statistics 


Figure  2.3  One-dimensional  N-grams 
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In  this  experiment,  binary  one-dimensional  sequences 
with  carefully-controlled  transition  probabilities, 
dependent  on  the  previous  four  points,  were  generated. 
Each  sequence  was  then  broken  up  into  shorter  512-pixel 
strings  which  were  stacked  to  form  the  two-dimensional 
512x512  array  which  served  as  the  texture  pattern  as  is 
illustrated  in  Figure  2.4.  Thus,  the  derived  statistics 
are  only  controlled  in  one  dimension  but  tne  final  texture 
is  two-dimensional. 


We  define  the  a  priori  probability  of  a  binary 
sequence  of  length  N  by  P  (V^, V2,  .  .  .  ,VN)  where  each  V^, 
i  =  1 , . . . , N  is  either  0  or  1.  In  our  experiment  this 
binary  sequence  is  determined  by  generation  parameters. 
This  set  of  parameters,  G  Q  (V  V  .  .  .  ,  ,  each  of  which 

represents  the  probability  of  generating  a  0  after  the 
contiguous  binary  sequence  V^, V^, . • • , V^,  defines  the 
Markov  process  used  to  generate  the  texture  pattern.  It 
follows  that  the  probability  of  generating  a  1  after  the 
sequence  V  ,V  , - 'VN  is  1-Go^Vl,V2' - 'VN)  •  That  is' 


vw* 


'V 


=  1-GX (V  fV2, 


V  )  (2.7) 
'  n' 


Illustration  of  this  commonly-used  texture  generation 
method  is  given  by  Purks  and  Richards  [3].  However,  it 
should  be  pointed  out  that  their  generation  parameters 
were  in  many  cases  constrained  to  provide  equal  N-gram 
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Figure  2.4  Forming  a  Two-Dimensional  Image  from 


statistics,  P (V^ , V2 , • . . , ) .  The  term  "N-gram"  refers  to 
the  Nth-order  joint  density  function  of  a  texture  and  was 
used  by  Purks  and  Richards  in  their  study  of  texture.  The 
term  MN-gram"is  frequently  used  in  this  thesis  as  a 
substitution  for  the  longer  terminology  "Nth-order  joint 
density  function."  A  texture  procedure,  more  general  than 
that  proposed  by  Purks  and  Richards,  is  detailed  here. 

These  generation  parameters  are  actually  a  special 
set  of  conditional  probabilities.  Only  this  specific 
group  of  conditional  probabilities  is  used  to  generate  the 
texture . 


2.3  Analytical  Analysis  of  the  Markov  Texture  Process 


By  examining  the  mathematics  of  the  Markov  process  we 
hope  to  generate  patterns  according  to  a  set  of  given 
probabilities  P (V^ , , . . . , )  which  may  be  named  the 
N-gram  statistics  of  a  specific  pattern.  We  must 
therefore  deal  with  the  relationships  that  exist  between 
these  N-gram  statistics  and  their  generation  parameters 
denoted  by  G (V^, • . • ,VN) .  Examining  these  relationships 
and  also  those  between  N-gram  statistics  of  different 
lengths  (that  is  the  relationships  between 


P (vl» V2' * ’ * »vn  )  and  p(Vi'V2, *  *  • ,VN2^  for  a11  N1  and  N2* 

leads  us  to  an  understanding  of  the  probabilistic  system 


involved  and  thereby  a  method  of  generating  desired 
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texture  patterns. 


In  generating  random  texture  sequences,  it  is  useful 

to  first  look  at  a  simple  analogy  to  the  process  from 

which  basic  concepts  and  conclusions  can  be  drawn.  We 

might  regard  this  generation  to  be  equivalent  to  the 

experiment  consisting  of  tossing  a  "smart"  coin  that  has  a 

finite  memory.  In  this  case,  Gq ' V2 ' ’ ' * ' Vn ^  might 

represent  the  probability  of  tossing  a  "heads"  given  the 

previous  sequence  of  N  tosses  was  V, ,V„,...,V  .  The 

12  N 

resulting  string  of  "l’"s  and  "0"'s  (0  is  the  random 

variable  denoting  heads,  1  denotes  tails)  recorded  from 

this  experiment  is  our  "texture."  We  realize  immediately 
that  the  texture  is  "determined"  by  this  set  of  generation 
parameters  G Q (V V . . . , VN) .  Using  the  concept  of 
conditional  probability  where  P(A/B)  is  the  probability  of 
A  given  B  we  notice  that 

G0(V1,V2"  '  *  ,VN}  =  P(0/V1,V2,...,VN)  .  (2.8) 

Perhaps  the  most  important  concept  derived  from  these 
generation  parameters  is  that  of  the  finite  memory  of  the 

system.  As  is  indicated  by  the  notation  Gq(Vj,  V2 . VN>' 

the  probability  of  generating  a  zero  depends  on  the  string 
of  binary  values  '  v2  '  •  •  •  »  vn  and  not  those  "preceding" 
•  It  is  thereby  suggested  that  our  system  has  an  N-gram 
memory  and  we  will  define  such  a  system  as 
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N-gram-dimensional .  For  example,  returning  to  our  coin 
tossing  experiment,  if  we  are  in  a  four-gram-dimensional 
system,  the  probability  of  tossing  a  head  depends  on  the 
four  previous  tosses  only  and  all  these  conditional 
probabilities  are  determined  by  the  sixteen  parameters 

G0(V1'V2'V3'V‘ 

With  these  concepts  in  mind  we  can  find  these 
generation  parameters  G  q(V  p  V2»  .  .  •  # VN)  given  the  desired 
probabilities  P (V^ , V2 , . . . , VN )  and  vice  versa.  The 
approach  taken  by  Purks  and  Richards  [3]  in  finding  these 
N-gram  statistics  is  based  on  sampling  the  generated 
textures.  This  may  be  seen  by  examining  the  entries  in 

t 

Table  I  of  Ref.  [3].  The  entries  correspond  to  the  number 
of  each  N-gram  counted  in  the  texture  generated  and  the 
accuracy  of  such  probabilities  depends  on  the  law  of  large 
numbers.  So  the  true  probabilities  P (V ^ , . . . , V^)  are 
only  approximated  by  the  output  textures  and  this 
approximation  is  poor  when  the  physical  size  of  the 
textures  is  small.  These  estimates  have  greater  variance 
when  the  true  probab i 1 i t i es  are  small.  Therefore  it  is 
desirable  to  compute  the  exact  probabilities  given  the 
generation  parameters  of  the  system. 

Before  proceeding  further  it  is  useful  to  prove  the 
identi ty 
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(2.9) 


P(V1'V2 . Vl1  = 


P(vrv2, . . .  »VN_1,0)+P(V1,V2, . . .  ,VN_1,1) 


Proof : 


P(v1,v2, . . . ,vN_lfi)  = 


P(VlfV2, . • •»VN_1) *GX (V1,V2, . . . ,VN_1)  (2.10) 


p(v1,v2, . . . #vN„1#o)  = 


P  ( ,  V2  /  •  •  *  » ^ )  *G0(VrV2,  .  .  .  / V^_  ^ ) 
P(V1,V2, . . . ,VN_1) * (1-G1 (VX,V2,. . . ,VN_1) ) 


therefore 


P(v1'v2* • • • /VN_1,0)+P(v1,v2, . . . ,vN_1,i)  = 

P<VV2 . VN-1}  * 

As  a  result  we  have  the  following  three  sets  of 
equal ities 

p(v1,v2, . . . ,vN_1,0)  = 

P(V1'V2» • • • » ^ ) *Gq  ^ vi ' V2 ' ’ *  * ' ^N  — 1 ^  (2.11) 

p(v1,v2, . . . ,vN_lfi)  = 

P  (V1,V2,  .  .  .  fVl).(l-G0  (V1,V2 . VN_X)  ) 


p(vrv2,...,vN_1)  = 


(2.12) 


P  (V^  ,V  2  t  •  •  • 


G0(V1'V2'  ‘  *  "VN-1)  =  p(V1,V2,  .  .  .  ,VN_1,0)/  (2.13) 

(P(vlfv2, . . . ,vN_1,0)+P(v1,v2, . . . ,VN_1,1) ) . 

Equation  (2.13)  results  from  Eqs.  (2.11)  and  (2.12)  and  is 
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essentially  a  statement  of  Bayes  theorem  for  our  problem. 

2.4  Texture  Statistics  From  Generation  Parameters 

Let  us  then  consider  the  problem  of  obtaining  the 
generation  parameters  (G's)  from  the  N-grams  (P's). 
Immediately  we  come  to  the  conclusion  that  this  is  a 
trivial  problem.  One  might  merely  use  Eq.  (2.13)  to 
deduce  the  generation  parameters.  However,  the  equations 
derived  so  far  do  not  indicate  the  extensive  relationships 
which  exist  between  the  N-grams. 

Equation  (2.13)  is  not  invertible  so  it  is  not  useful 

in  obtaining  the  N-gram  statistics,  P(V^,...,VN)  of  a 

sequence  given  the  generation  parameters.  As  was  stated 

above,  once  the  generation  parameters  are  defined,  a 

texture  may  be  generated  using  those  parameters  and  the 

N-gram  statistics  are  determined.  We  also  know  that  once 

a  complete  set  of  N-gram  statistics  P(V..,...,V  )  are 

N1 

defined  for  some  N  ,  the  N-gram  statistics  P  (V  ,...,V  ) 

1  y  1  n2 

may  be  resolved  using  Eq .  (2.12)  for  all  N2<Ni*  Given  the 

generation  parameters  of  a  system,  G ^ (V^ , V2 , . . . , ) ,  we 

can  analytically  determine  the  N-gram  statistics, 

P  (Vi ,  ^2  '  *  *  *  » VN )  resulting  texture. 

The  solution  to  this  problem  of  finding  N-gram 
statistics  given  generation  parameters  may  be  found  by 


24 


considering  the  generation  procedure  as  a  discrete  state 

Markov  process.  This  approach  is  readily  seen  wnen 

considering  the  generation  parameters  G- (V,  ,  V_  ,  .  .  .  ,  V  )  as 

u  1  2  N 

transition  probabilities.  If  we  consider  a 

two-dimensional  system  with  P  (0 , 0  )  ,  P  (0  ,  1 )  ,  P  (1 , 0  )  and 

P(l,l)  and  generation  parameters  G(0,0),  G(0,1),  G(1,0) 

and  G(l,l)  we  may  define  our  system  as  composed  of  four 

possible  states  (0,0),  (0,1),  (1,0)  and  (1,1).  If  the 

system  is  in  state  i  at  the  Kth  observation  and  in  state  j 

at  the  (K+l)th  observation  then  we  say  that  the  system  has 

made  a  transition  from  state  i  to  state  j  at  the  Kth  stage 

of  the  generation  process.  In  our  example  an  observation 

is  taken  at  each  generation  of  a  single  new  binary  value 

and  the  state  is  determined  by  the  values  of  the  last  two 

binary  numbers  generated.  As  an  example,  consider  the 

sequence  0,1, 1,0,0.  We  might  consider  the  system  to  be  in 

the  (0,1)  state  at  the  start  which  may  represent  the  Kth 

stage  of  our  generation  process  then  a  transition  is  made 

to  the  (1,1)  state  at  the  (K+l)  stage.  These  transition 

probabilities  are  determined  by  the  generation  parameters 

of  the  system.  We  also  note  that  our  N-dimensional  system 
N 

has  2  possible  states.  As  the  transitions  from  each  of 
these  2N  possible  states  to  each  of  the  2N  possible  states 
is  fixed  by  our  generation  parameters  we  may  form  a 
transition  matrix  T  whose  elements  t(i,j)  represent  the 
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probability  of  a  transition  from  the  ith  state  to  the  j th 
state.  If  T  is  the  transition  matrix  of  a  regular  Markov 
chain,  then  there  is  a  unique  probability  vector  p  which 
has  positive  coordinates  and  satisfies 

TTp  =  p  (2.14) 

This  same  vector  p  may  be  computed  by  taking  any  row  of 
the  matrix 

Tq  (2.15) 

as  q  approaches  infinity  [17].  The  vector  p  represents 
the  vector  of  steady  state  probab i 1 i t i es .  In  our  case  it 
contains  the  desired  probabilities  P (V^ , V2 , . . . , VN )  . 

It  is  important  to  realize  that  this  theorem  holds 
for  regular  Markov  processes.  If  there  is  an  integer  q 
such  that  every  element  of  the  matrix  in  Eq .  (2.15)  is 
stricly  positive  then  the  process  is  regular.  Some 
processes  are  not  regular  such  as  absorbing  Markov  chains 
[17].  When  any  element  of  the  transition  matrix  is  equal 
to  one  along  the  diagonal,  the  process  is  said  to  be 
absorbing  given  that  the  system  may  begin  in  any  state. 
This  could  happen  if  Gg(0,0)  =  1  for  example  (a  series  of 
0's  would  be  generated  in  this  case).  For  the  purposes  of 
our  discussion  we  will  assume  that 


26 


for  all  V-^.  This  is  a  sufficient  but  not  necessary 
condition  for  the  process  to  be  regular. 


Appl yi ng 

these 

concepts  to 

a  two- 

d  imens ional 

system 

we  obtain  the 

transition  matrix 

Final 

State 

- 

00 

01 

10 

11 

00 

o 

o 

o 

o 

1-Gq (0,0) 

0 

0 

Initial 

State 

01 

0 

0 

Gq(0,1) 

1-Gq (0,1) 

(2.17) 

10 

O 

r— H 

o 

o 

1-G0(1,0) 

0 

0 

.11 

0 

0 

G0(l,l) 

1-G0(1,1)_ 

The  first  row  contains  the  transition  probabilities 
from  state  (0,0)  to  states  (0,0),  (0,1),  (1,0)  and  (1,1) 
in  that  order.  The  following  set  of  equations  results 
when  Eq.  (2.14)  and  Eq .  (2.17)  are  combined: 


GQ  (0 , 0 ) -1 

0 

Gq  (1,0) 

0 

P(0,0) 

0 

1-Gq (0,0) 

-1 

1-G0(1,0) 

0 

X 

P(0,1) 

0 

0 

GQ  ( 0 , 1 ) 

-1 

G0(l,l) 

P(1,0) 

0 

0 

1-Gq (0,1) 

0 

-G0d,l) 

P(U) 

_0_ 

As  the  above  system  is  singular ,  we  may  form  an 


equivalent  non-singular  set  by  replacing  any  equation  with 


P(0,0)+P(0,1)+P(1,0)+P(1,1)  =  1  (2.19) 

by  using  the  fact  that  p  is  a  probability  vector.  Solving 
this  system  gives  the  desired  two-gram  statistics 

p(v1,v2) . 

Examining  these  generation  parameters  further  we  find 
that  the  same  N-gram  statistics  may  be  generated  by 
generation  parameters  of  a  different  dimension.  For 
example,  consider  the  generation  parameters  in  Table  2.1. 


TABLE  2.1.  EQUIVALENT  GENERATION  PARAMETERS 


Set  1 

Set 

_3 

Gq  (0,0)  =  0.05 

o 

o 

o 

o 

=  0.05 

Gq  (1,0,0)  =  0.05 

Gq  (0,1)  =  0.07 

Gq (0,0,1) 

=  0.07 

Gq(1,0,1)  =  0.07 

Gq(1,0)  =  0.92 

Gq(0,1,0) 

=  0.92 

G  (1,1,0)  =  0.92 

Gq  (1,1)  =  0.75 

Gq(0,1,1) 

=  0.75 

Gq(1,1,1)  =  0.75 

Notice  that  the  probability  of  generating  a  zero 

following  a  vi»v2'v3  ^oes  not  depend  on  V^.  The  values,  | 

| 

Go(Vi,V2,V3)  of  the  second  set  indicate  that  the  system  is-  } 

memoryless  beyond  two  previous  generation  steps.  We  may 

write 
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G0(V1,V2'V3)  =  P(0/VlfV2,V3) 


(2.20) 


p(o/v2,v3)  =  g0(v2,v3) 

for  all  V2  and 

It  follows  that,  according  to  Eq.  (2.11),  the 
P  (V^  ,  V  ,  .  .  .  ,V  )  are  also  determined  in  our  example  for  N>2 
given  P(V^,V2)  and  gq(V]_»V2^*  Thus  we  have 


p(v1,v2,v3)  =  p(v1,v2)*gv  (v1#v2) 


(2.21) 


P(V1'V2,V3,V4)  =  P(V1'V2'V3) *GV  (V1'V2'V3} 

=  P(V1,V2,V  ).G v  (V  ,V  > 

4 


(2.22) 


We  conclude  that  given  any  set  of  generation  parameters 

G0(Vi,V2 . VN)  we  may  form  a  set  of  generation 

parameters  Gg  (V-]_ ,  V2  ,  .  .  .  ,  VN )  ,  M  greater  than  or  equal  to  N, 
according  to  the  rule 


G0  {V1,V2'  '  •  ’  ,VM-N'VM-N+1'  *  *  •  'V  “ 

go(Vn+i . V 


(2.23) 


that  generate  an  equivalent  set  of  N-gram  statistics  and 
therefore  equivalent  textures. 


When  desiring  to  generate  textures  according  to  a 
given  set  of  N-gram  statistics,  P  (V^ , V2 , . . . , V  )  we  must 
realize  the  set  of  constraints  imposed  on  the  set.  For 
example , 
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(2.24) 


r 


'V 


1 


for  all  N.  Returning  to  the  set  of  equations  used  to 
determine  the  2-gram  statistics  in  matrix  form  we  realize 
that,  by  adding  the  first  two  rows  and  last  two  rows  of 
Eq.  (2.18),  P(0,1)  =  P(1,0).  In  fact,  by  considering  the 
set  of  equations  arising  from  the  set  of  equations  derives 
from  the  generation  systems  of  higher  dimensions  we  find 


P(V2 


,v2,v2,v2, . 
p(v3,v2,v2 


'•'W  = 


,v2, . . . ,v2,vx) 


(2.25) 


This  implies  that  many  constraints 
N-gram  statistics.  For  example,  in 
system  containing  P(V^,V2,V3) 

P(0,0,1)  =  P  (1 , 0 , 0) 
P(0,1,1)  =  P ( 1 , 1 , 0) 


are  present  on  the 
the  3-gram-dimensional 


(2. 26) 


but  also  by  Eq.  (2.12)  and  the  fact  that  P(0,1)  =  P(1,0), 


P(0,1,0)+P(0,1,1)  =  P(1,0,0)+P(1,0,1) 


(2.27) 


and 


P(0,0,1)+P (1,0,1)  =  P(0,1,0)+P (1,1,0) 


By  definition  we  also  know  that 
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P(0,v1,v2,- . . ,vn_1)*g0(0/v1,v2, . . . ,vN_1) 


+  P(1,V1,V2,  .  .  .  /VN_1)  *GQ(1/V1,V2,  .  .  .  /V^)  (2.28) 

=  p(v1,v2, . . . /VN_1,0) 


and 


p(0,v1,v2, 

+P(1/V1,V 
=  P(VlfV 


2 

2 


'Vi)*[1-G0(0'V1'V2'‘ 

••'VN-1)*[1_G0(1'V1'V2 


'Vl)] 


•  ’ ,VN-1) 1 


(2.29) 


Combining  Eqs.  (2.28)  and  (2.29)  with  Eq .  (2.12) 

P(0,V1'V2' * ’ * ,VN-1)+P(1,V1'V2' ‘ ' * 'VN-1)  = 

p(vlf . . . ,vN_x) 

Equation  (2.30)  holds  for  all  N>2. 


(2.30) 


Four  constraining  equations  exist  when  Eqs.  (2.24), 

(2.25)  and  (2.30)  are  reduced  to  orthogonal  form  for  this 

3-gram-dimensional  system.  This  implies  that  we  have  four 

degrees  of  freedom  when  choosing  the  eight  values 

P^1,V2'V3^  along  with  the  obvious  constraint  that 

P (V^ , ^2  > v3 ) >0 •  This  is  precisely  the  number  of  degrees  of 

freedom  in  the  set  of  gq^V1'V2^  which  have  only  the  range 

constraint  of  Eq.  (2.16).  For  higher  N-grain-dimensional 

N-l 

systems ,  P  (V-^ ,  . .  .  'V  always  has  2  constraints  on  it 

from  Eqs.  (2.12),  (2.24),  (2.25)  and  (2.30).  Thus  we  see 
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that  Eq .  (2.13)  holds  in  a  degrees  of  freedom  sense. 

In  conclusion,  a  method  of  determining  N-gram 
statistics  from  generation  parameters  using  the  concept  of 
a  Markov  chain  was  developed  and  a  set  of  linear  equations 
describing  the  constraints  on  the  N-grams  was  presented. 
Using  this  approach  to  generating  texture  patterns  a  large 
variety  of  textures  may  be  easily  generated  and  examined 
using  a  minimum  amount  of  effort. 

2.5  Texture  Moments 

The  above  equalities  and  inequalities  provide  a  full 

understanding  of  the  texture  generation  system  in 

probabilistic  terms.  Still  further  conclusions  can  be 

derived  from  them.  From  the  above  we  see  more  clearly 

that  the  generation  parameters  G  (Vj_,  V2,  . . .  ,VN)  determine 

the  texture  completely  and  thus  define  the  N-gram 

statistics  P(V^,...,VM)  for  all  M.  Also  for  a  given  set 

of  P (V . , . . . , V  )  there  can  exist  an  infinite  number  of 
i  M 

generation  parameters,  G (V^, V2, . . . ,V  ) ,  which  would 
generate  many  textures  with  such  statistics  if  N>M-1. 
Provided  the  constraints  on  the  statistics  P(V^,...,V  ) 
are  met  just  one  texture  could  be  generated  if  N  =  M-l  or 
perhaps  none  at  all  if  N<M-1.  Thus  textures  with  equal 
first,  second,  third  and  fourth  nearest-neighbor 
probabilities  can  be  generated. 
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Parameters  thought  to  be  useful  in  texture 
discrimination  may  also  be  easily  developed.  For  example, 
joint  moments  about  the  mean  defined  as 


1C  T  T  T 

E[(x1-u1)  1(x2-u2)  2(x3-y3)  3...(xk-,k)  k]  (2.31) 


y* 

r  .  i! 


where  ^  r^  is  the  order  of  moment  [18].  The  rth  moment 
of  x^  is  defined  as 


E  (x^)  x^f  (x1#  .  .  .  ,xk)  (2.32) 

X1  X2  Xk 

where  f  is  the  joint  probability  distribution  of  the  x^. 
From  our  binary  textures  we  could  define  the  following 
parameters : 

u  =  E{Xo}  =Xlxif(xi)  =  SXi  P(xi}  = 

0-P(0)+l-P(l)  =  P  (1) 
a2  =  E{  (Xq~u  )  2  }  =^^(x-y)  2f  (x)  = 

(O-P(l)  )2P(0)  +  (1-P(1)  )2P(1) 

=  P(l)-P(l)2 

E  {  (xQ-y )  3 }  =^(x-y)3f(x)  =  (2.33) 

(O-P(l)  )3P(0)  +  (1-P(1)  )  3P(1) 

-2P(1) 3-3P (1) 2+P (1) 

E{  (xQ-p  )  (Xj-y  )  }  _  P  (1 , 1)  -P  (1)  2 

a  —  2  —  ■  5 

a  P(l)-P(l) 


33 


1 


P ( 1 , 1 , 1) -3P ( 1 , 1)  *P  ( 1 ) +4P (1 ) 3 
(P(l)-P(l)2)3/2 

where  P(l),  P(l,l),  P(l,l,l)  represent  nearest-neighbor 
(N-gram)  statistics  although  this  can  be  changed  to 
include  non-nearest-neighbor  statistics  thus  creating  new 
texture  parameters.  The  above  parameters  are  useful  in 
d iscr im inat ion  therefore  only  when  textures  differ  in 
their  (3-gram)  or  shorter  statistics. 

2.6  Constraining  Second-Order  Statistics 

We  describe  now  a  method  which  allows 
non-nearest-neighbor  statistics  to  be  controlled  using  the 
relationships  developed  in  the  Sections  2.3  and  2.4. 
Because  second-order  probabilities  are  of  interest  we 
investigate  the  conditions  required  to  assure  equality  of 
second-order  statistics  for  non-nearest-neighbor 
statistics.  If  we  denote  G  g  (V  V  2»  .  .  •  ,  VM)  as  our 
generation  parameters  and  P  (V  V  2,  •  •  •  ,  VN)  to  be  their 
associated  N-gram  statistics  then  fo'r  the  second-order 
statistics  of  one  texture  to  be  equal  to  another  for  the 
(N-l)st  neighbor  distance, 


E{ (x0-y) (x1~u) (x2~u) } 

3 

a 


3 


(2.34) 


X  •••X)  Pa<VV2 . VVM+1 . V 


v2  VN-1  (2-34) 

E-  ’  '  EPb<Vl'V2"-- ,VM+1'  •  •  •  ,VN} 

'  V2  VN-1 

for  V-^  ,  VN  e{0 , 1 }  where  Pa  ,  represents  the  N-gram 

statistics  for  the  first  and  second  texture  respectively. 
It  can  also  be  shown  that 


E  •  E  p(vi'v2 . vm . v  = 


V2  VN-1 


E  **  *  E  p(vi,v2'  *  •  • ' 


v 


(2.35) 


V2  VN-1 


V, 


where 


V  E.,IV("1)  k  •  G0<Vk-M'Vk-M+l . Vl11 

k=M+l 


Vj  Vj=0 

Recall  that  the  N-gram  statistics,  P,  are  a  function 
of  the  generating  parameters  G.  If  the  second-order 
statistics  are  to  be  equal  for  two  textures  regardless  of 
neighbor  distance  then  Eq .  (2.35)  must  hold  for  all  N.  It 
was  previously  believed  that  combining  Eq .  (2.34)  and 
Eq.  (2.35)  yielded  a  non- redundant  non-linear  set  of 
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distributions  cannot  exist.  He  then  proceeded  to  work 
with  textures  of  four  grey  levels. 


In  the  next  section,  we  will  show  that  binary 
textures  with  equal  second-order  statistics  for  all 
distances  can  be  generated  using  a  carefully-chosen  set  of 
generation  parameters  and  that  these  texture  contadict  the 
Julesz  conjecture. 

2.7  Experimental  Results 

We  may  define  second-order  statistics  for 
non-nearest-neighbors  as 


(2.36) 


p<w  -Z  •••X)  Z  •••Ep(vv2 . vj . 


V 


V. 


Vj-1  Vl  VN 


Purks  and  Richards  [3]  attempted  to  demonstrate  that 
textures  equal  in  second-order  distribution  could  be 
generated  with  visually  detectable  differences.  However, 
they  merely  held  second-  or  third-order  statistics  equal 
between  the  two  textures  for  small  j,  while  varying 
second-  and  third-order  statistics  for  longer  j,  that  is, 
they  allowed  the  second-order  statistics  for 
non-nearest-neighbors  to  be  unequal  over  some  distance. 
Examples  of  this  type  of  manipulation  are  shown  in 
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Figs.  2.5-2.10.  The  corresponding  N-grams  are  shown  in 
Table  2.2  and  Table  2.3.  The  first  column  of  each 

subtable  contains  the  analytic  N-grams  for  texture  (a), 
top,  and  texture  (b),  bottom,  which  are  found  by  solving 
Eq.  (2.14)  and  Eq .  (2.24)  given  the  generation  parameters 
in  the  third  column.  The  middle  column  contains  the 
actual  N-gram  statistics  as  measured  from  the  generated 
textures  themselves.  The  input,  analytically-solved 

parameters  will  not  be  equal  to  the  output  statistics  as 
the  generation  process  is  random.  In  other  words,  the 
output  statistics  are  based  on  measurements  taken  on  a 
sample  from  a  population  of  textures  having  the 

characteristics  defined  by  the  input  parameters. 

Figures  2.5  and  2.6  show  pairs  of  statistics  having 

equal  first-  and  second-order  nearest-neighbor  statistics. 

That  is,  P  (V  , V  )  =  P,  (V,,V„).  There  is  also  an  internal 
a  1  2  b  1  2 

equality  for  these  textures  which  may  be  expressed  as 

pa (Vf , v2 )  =  Pb^vl'v2)  =  1/4  for  a11  V1  and  v 2  that  are 

nearest-neighbors.  Visual  differences  are  apparent 

between  the  pairs.  Figures  2.7  and  2.8  have  texture  pairs 
which  have  equal  third-order,  nearest-neighbor  statistics 
both  within  and  between  the  pairs,  that  is, 

Pa(V1,V2,V3)  =  Pb(VlfV2,V3)  =  0.125.  The  bottom  texture 
in  Fig.  2.7  is  a  coin  flip,  purely  random  texture  with  the 
probability  of  both  black  and  white  equal.  It  is 
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Table  2.3  N-GRAMS  AMD  GEHERATIOM  PARAMETERS 
FOR  FIGURES  2.9  -  2.12 
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interesting  to  note  that  this  texture  is  visually  quite 
similar  to  the  bottom  texture  of  Fig.  2.6,  yet  the 
generation  parameters  and  N-grams  differ  significantly. 
Thus  two  textures  with  differing  N-gram  statistics  may  be 
generated  which  are  visually  similar.  Figure  2.9  shows  a 
pair  of  textures  which  .  have  equal  fourth-order, 
nearest-neighbor  statistics  both  within  and  between  the 
pairs.  That  is, 


Pa(Vl'V2'V3'V4)  =  Pb(Vl'V2'V3'V4)  =  °'0625  •  (2'37) 


The  two  are  easily  discriminated.  Figure  2.10  shows  two 
textures  with  between-equal ,  fourth-order  nearest-neighbor 
statistics.  The  N-grams  are  not  equal  within  however. 


Figures  2.11  and  2.12  have  texture  pairs  which  are 
similar  in  many  ways  to  the  pair  in  Fig.  2.6.  Both  are 
counter  examples  to  Julesz's  conjecture  that  the  eye  is 
sensitive  to  only  first-order  and  second-order  probability 
distributions.  Each  has  a  texture  pair  where  second-order 
statistics  are  equal  for  all  distances  between  the  texture 
pair.  Precisely  stated,  pa(vi»vj)  =  pb^vl,vj^  f°r  a33  3* 
These  textures  are  generated  using  a  set  of  restrictions 
discussed  in  Appendix  A. 


The  textures  in  these  figures  are  pseudo-randomly 
generated.  Each  texture  was  tested  using  a 
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goodness-of-f it  procedure  to  insure  that  textures  with  the 
desired  N-grams  were  generated.  This  was  done  for  1-grams 
to  10-grams.  The  procedure  is  outlined  in  Appendix  B. 
The  goodness-of- f it  depends  on  the  pseudo- random  number 
generator  used,  as  a  poor  generator  can  be  guaranteed  to 
yield  poor  results  in  this  type  of  experiment.  The 
pseudo- random  number  generator  used  is  detailed  in 
Appendix  C. 

2.8  Conclusions 

The  one-dimensional  binary  patterns  generated  give 
rise  to  some  basic  concepts  concerning  textures  and  their 
discrimination.  First  of  all  they  indicate  the  use  of 
moments  and  similar  statistics  is  not  optimal  at  least  in 
the  nearest-neighbor  sense  as  many  textures  have  equal 
moments  but  are  visually  quite  different.  However,  it 
should  be  pointed  our  that  this  may  only  be  characteristic 
of  some  artificial  textures  and  that  moments  could  serve 
as  good  discrimination  parameters  in  many  real-world 
applications.  Secondly,  the  results  indicate  a  close 
relationship  between  second-order  non-nearest-neighbor 
statistics  and  human  discrimination.  The  counter-examples 
to  the  Julesz  conjecture  indicate  that  N-grams  and 
higher-order  statistics  may  be  valuable  in  identifying  and 
discriminating  some  textures.  Use  of  these  texture 
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measures  depends  on  factors  such  as  discrimination 
accuracy  desired,  cost  factors  for  statistics  measure  and 
the  nature  of  the  textures  involved. 


In  later  chapters,  investigation  of  two-dimensional 
textures  will  be  pursued  as  these  correspond  more  closely 
to  natural  scene  textures.  The  texture  generation  task 
will  be  approached  from  a  simulation  rather  than  a  pure 
synthesis  point  of  view.  The  complexity  of  controlling 
two-dimensional  statistics  in  a  synthesis  process  and 
attempts  to  study  their  effects  on  human  discrimination 
and  interpretation  of  textures  is  a  problem  which  requires 
careful  analysis  and  is  beyond  the  scope  of  this  work.  It 
is  very  possible  that  no  applicable,  clear-cut  results 
could  be  obtained  from  such  a  study.  However,  by 
simulating  textures,  models  and  processes  for  generating 
similar-looking  textures  are  derived  which  may  be  useful 
in  other  applications  and  some  simple  statements 
concerning  human  discrimination  of  textures  can  also  be 
made . 
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CHAPTER  3 


TWO-DIMENSIONAL  BINARY  TEXTURE  MODELS 

3.1  Introduction 

In  this  chapter,  the  concepts  used  to  generate 
one-dimensional  binary  textures  are  extended  to  the 
two-dimensional  case.  This  model  is  then  used  to  simulate 
natural  binary  textures.  A  method  for  choosing  the  pixels 
in  a  non-cont ig uous  generation  kernel  based  on  a  linear 
model  is  described.  This  is  an  important  concept  in  much 
of  the  work  presented  in  this  thesis.  The  method  for 
collecting  N-gram  statistics  is  discussed  and  practical 
problems  arising  in  this  process  are  investigated. 
Results  of  natural  texture  simulation  using  this  method 
are  presented.  Finally,  the  linear  model  which  was  used 
to  determine  the  kernel  pixels  of  greatest  value  in  the 
N-gram  generation  process  is  used  to  generate  binary 
simulations.  This  will  lead  us  to  the  application  of  the 
linear  model  to  continuous-tone  textures  in  Chapter  5. 

3.2  The  2-D  Binary  Markov  Model 

In  the  investigation  of  natural  phenomenon  once  a 
researcher  collects  enough  data  he  tries  to  imagine  a 
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process  :ich  accounts  for  the  results.  The  construction 
and  development  of  a  mathematical  model  is  often  the  best 
way  to  do  this.  done.  In  some  cases,  the  model  may  be 
extraordinarily  complex,  in  others,  exceedingly  simple. 
In  most  cases  model  acceptance  cannot  be  based  on  "truth" 
as  the  true  generating  phenomenon  is  too  complex  or  simply 
unknown  and  so  it  is  based  upon  model  usefulness  and  "how 
well  it  works".  Such  "working"  models  for  texture  are 
presented  in  this  chapter  because  they  have  the  ability  to 
simulate  some  natural  textures. 

If  one  can  synthesize  and  simulate  natural  textures 
adequately  by  using  some  proposed  model  the  criterion  of 
usefulness  and  workability  for  that  model  is  met.  A 
researcher  may  then  also  apply  the  model  to  problems  of 
texture  identification  and  discrimination  with 
justification.  With  any  set  of  texture  measures  required 
for  simulation,  the  information  content  of  the  measures  is 
viaually  indicated  by  the  quality  of  the  simulation  and 
not  merely  by  the  percent  of  correct  versus  incorrect 
classifications  when  those  textures  are  applied  to  the 
d isc r im ina t ion  problem.  By  adding  features  important  to 
texture  simulation  better  methods  of  texture 
discrimination  and  identification  can  possibly  be  found. 

Many  early  texture  studies  involved  the  use  of  binary 
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textures  generated  by  one-dimensional  Markov  processes. 


Such  work  was 
one-dimensional 
generated  line 
pa  rameter s 


presented  in  Chapter  2. 
models  a  large  vector  of 
by  line  using  a  set  of 


In  these 
pixels  was 
generation 


where 


V  (V1'V2 . V 

N+l 


GV  {V1,V2'  -  ‘  *  ,VN>  =  P(VN+1/V1'V2'-"'VN)  (3,1) 
N+l 

and  P(A/B)  represents  the  probability  of  A  given  B.  In 
the  above  notation  each  \A  represents  a  generated  pixel 
which  has  value  0  (black)  or  1  (white).  Each  pixel  value, 
then,  depends  on  the  N  pixels  previous  to  it.  A 
two-dimensional  texture  image  is  then  formed  by  breaking 
up  the  large  vector  of  pixels  into  shorter  strings  and 
stacking  them  one  on  top  of  the  other  (see  Fig.  2.4). 
This  procedure  for  large  images  nearly  insured  image  row 
independence  (unless  N  was  large)  thereby  creating  only 
horizontally  oriented  textures  totally  unsuitable  for 
simulating  natural  two-dimensional  textures. 


By  allowing  N  to  increase  exceeding  the  short  string 
line  length,  two-dimensional  (vertical  and  horizontal) 
dependence  may  be  induced  into  the  generating  process.  A 
pixel  value  then  depends  not  only  on  the  pixels  previous 


to  it  on  the  same  line  but  also  on  the  pixels  above  it 

(see  Fig.  3.1(b)).  Thus,  textures  could  be  generated  as  a 

time  sequence  in  television  raster  scan  fashion.  In 

theory,  texture  dependence  could  be  extended  _ad  infinitum, 

however  practical  considerations  concerning  the  actual 
•  N 

generation  process  show  us  that  2  generation  parameters 
must  be  accounted  for.  As  a  possible  solution  to  the 
storage  problem  we  can  choose  to  ignore,  all  but  N  of  the 
previous  pixels  in  our  generation  process  and  we  can  allow 
the  pattern  of  the  ' s  to  become  flexible.  This  idea 
will  be  discussed  in  later  sections  of  this  chapter. 

Throughout  the  remainder  of  this  thesis,  the  set  of 
pixels,  VVs,  on  which  the  next  pixel,  VN+1,  depends  will 
be  referred  to  as  the  "kernel"  of  the  synthesis  process. 
The  pixel  VN+1  will  be  referred  to  as  the  "eye"  of  the 
kernel . 

In  order  to  estimate  P  (V^^,  V2,  . .  ,VN)  for  a  fixed 
pattern  V^, V^, . . . , V^,  all  M  substrings  (samples)  of  length 
N  are  taken  from  a  parent  substring  of  length  M+N-l  and 
the  number  of  occurrences  of  the  specific  pattern 
V  ,V  , ...,V  are  counted,  then  divided  by  M.  This  is 
equivalent  to  estimating  the  probability  density  function 
of  a  random  variable  by  the  histogram  of  a  set  of  samples. 
Ignoring  boundary  conditions,  linear  unbiased  estimates 
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(P's)  of  the  P's  may  be  defined  as 


P(V1,V2 


'V 


1 

M 


N 

n  6 (I (k+j ) -v, ) 

k=l  K 


(3.2) 


where 
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1 


if  Vj  *  Vk 


if 


V  . 
D 


V 


k 


and  I ( i )  represents  the  ith  element  of  the  one-dimensional 
texture  string  from  which  the  parameters  are  to  be 
estimated.  Equation  (3.2)  assumes  that  the  V.  are 

l 

contiguously  located  in  order  along  a  line. 


This  idea  of  estimating  N-grams,  P(V_  ,  V  .  .  .  .  ,V  )  from 

12  N 

a  sample  parent  texture  may  be  extended  to  the 
two-dimensional  case.  A  histogram  of  occurrences  of  each 
pattern  of  (V^ , , . . . , )  is  made  by  passing  the 
two-dimensional  kernel  in  Fig.  3.1(b)  over  the 
two-dimensional  sample  parent  image.  The  tally  is  then 
divided  by  the  total  number  of  sample  patterns  observed  to 
obtain  P(V^f...,V  ).  As  was  stated  earlier, 

two-dimensional  synthesis  is  merely  an  extension  of  the 
one-dimensional  case  ignoring  boundary  conditions  of  the 
two-dimensional  image. 


Although  it  was  not  explicitly  stated  in  Chapter  2, 
the  generation  parameters  of  a  texture  may  be  estimated 
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These 


for  any  given  set  of  N  V^'s  from  a  parent  texture, 
statistics  have  the  property 


E  [G 


V  (vi » v2 ' 

N+l 


,V  ) ]  =  G  (V  ,V  , . . . ,VM) 
N  VN+1  1  2  N 


where 


V  ,  (V1-V2 . V  ■  P(VV2 . VN-Vl)/ 

N+l 

(P(v1,v2 . VN,0)+P(V1,V2 . V  })  . 


(3.3) 


3.3  Seeding  the  Generation  Process 

In  the  one-dimensional  texture  pattern  generation  of 
Chapter  2,  a  long  binary  vector  of  pixels  is  broken  up 
into  short  strings  which  are  stacked  on  top  of  each  other 
(see  Fig.  2.4).  When  synthesizing  this  one-dimensional 
string,  a  seed  of  four  binary  values  is  actually  required 
to  begin  the  process  for  a  4-gram-dimensional  system. 
When  the  Markov  process  is  regular  and  non-absorbing  with 
none  of  the  generation  parameters  equal  to  one  or  zero  the 
process  rapidly  reaches  a  steady  state  independent  of  the 
seed . 


A  similar  thing  happens  in  the  two-dimensional  case 
but  here  the  seed  is  larger  and  forms  a  two-dimensional 
frame  around  the  synthesized  image  texture  as  shown  in 
Fig.  3.2.  With  a  two-dimensional  texture  generation 
kernel  such  as  that  shown  in  Fig.  3.1(b),  the  synthesized 
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image  texture  may  be  generated  without  using  the  bottom  of 
the  frame  so  the  next  pixel  at  each  generation  step 
depends  only  on  pixels  above  it. 

The  seeding  process  may  be  handled  in  a  variety  of 
ways.  The  simplest  approach  would  be  to  randomly  generate 
the  seed  once  for  the  whole  image.  In  this  case  the  pixel 
values  in  the  seed  frame  of  Fig.  3.2  remain  the  same 
throughout  the  generation  process.  A  second  approach 
would  require  the  random  generation  of  each  pixel,  V.,  in 
the  generation  kernel  that  fell  outside  the  synthesized 
image  texture  region  at  each  step  during  the  generation 
process.  Using  this  method,  the  pixel  values  in  the  seed 
frame  change  at  each  pixel  generation.  A  third  method 
would  involve  wrapping  the  image  around  such  that  the  left 
edge  of  the  synthesized  image  texture  joins  with  the  right 
side  as  in  Fig.  3.3.  In  this  case  a  random  seed  is 
required  only  to  begin  the  process  of  the  top  of  an  image. 
A  final  method  involves  the  use  of  another  texture, 
usually  a  previously  generated  texture  or  the  parent 
texture,  as  part  of  the  seed,  rather  than  noise.  This 
method  reduces  the  noise  which  otherwise  occurs  around  the 
edges  of  the  synthesized  image  texture. 

Regardless  of  the  seeding  process,  all  texture 
synthesis  methods  developed  in  this  thesis  normally 
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converge  to  a  steady  state  within  5  to  20  pixels  of  the 
border  of  the  image.  This  was  confirmed  by  repeated 
studies  of  convergence  effects  on  texture  simulations.  In 
most  cases,  this  narrow  region  is  not  noticeable  and  is 
included  as  a  part  of  the  result.  In  some  critical 
applications  these  edges  could  be  thrown  away. 

3.4  Kernel  Selection  Using  The  Linear  Model 

We  will  refer  to  the  V^'s  on  which  the  next  pixel, 
V  depends  as  the  kernel  of  the  generation  process. 

Geometr ical ly  speaking,  the  V^'s  form  a  kernel  "shape"  or 
"pattern"  which  may  or  may  not  be  spatially  contiguous. 
For  example,  in  Fig.  3.4  a  generating  kernel  shape  is 
shown  where  the  V5  pixel  directly  depends  on  only  some  of 
the  pixels  in  its  surrounding  neighborhood.  In  this  case, 
V,_  may  be  generated  based  on  the  values  of  pixels 
Vf,  V  ,  V3  and  V4  but  is  directly  dependent  on  no  other 
pixels  in  the  neighborhood.  This  does  not  imply  that 
is  not  related  or  correlated  with  its  other  neighbors.  In 
fact,  the  relationships  between  V^,  V 2»  V^,  V^,  and 
will  determine  other  interrelationships. 

A  non-contiguous  neighborhood  of  V^'s  is  used  as  it 
allows  a  more  parsimonious  model  for  texture  generation  to 
be  chosen.  An  analogy  is  in  simple  linear  regression  (as 
defined  by  Draper[20])  where  independent  variables  which 
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do  not  contribute  to  the  prediction  or  estimation  of  the 


dependent  variable  are  dropped.  In  texture  generation 

this  allows  the  model  to  be  estimated  by  fewer  parameters 

and  makes  the  generation-synthesis  process  more  efficient 

by  reducing  the  number  of  computations  required.  When 

generating  textures  based  on  N-grams,  reducing  N  reduces 

N 

the  amount  of  storage  required  for  2  generation 

parameters.  By  allowing  the  kernel  of  ' s  on  which  VN+^ 

depends  to  be  non-contiguous ,  the  range  of  dependence  in  a 
distance  sense  is  increased  over  that  which  would  be 
allowed  with  a  contiguous  kernel  containing  the  same 
number  of  V^’s.  This  is  very  important  to  obtain  the 
larger  structure  apparent  in  many  textures.  Reducing  the 
number  of  pixels  in  the  model  also  relieves  us  from  the 
complex  numerical  problems  of  inverting  matrices  of 
unwieldy  size,  a  necessary  step  in  linear  model  parameter 
estimation  discussed  later  in  this  chapter.  We  would,  for 
example,  not  expect  our  VN+^  pixel  to  depend  on  a  pixel 
where  the  spatial  separation  between  VN+^  and  is  large. 
If  that  distance  is  small,  however,  we  would  expect  a 
large  dependence. 

The  method  for  choosing  the  proper  independent 
variables  (V^'s)  to  be  included  in  the  generation  process 
requires  special  attention.  We  wish  to  choose  the  best 
subset  of  N  variables  from  a  larger  finite  neighborhood  of 


T  variables,  where  N<T .  Evaluating  such  subsets  and  their 

corresponding  models  requires  a  criterion.  Texture 

results  for  each  possible  model  could  be  visually  examined 

and  compared  and  the  V^’s  of  the  model  corresponding  to 

the  visually  most  pleasing  result  could  be  chosen. 

T 

However,  (N )  model  evaluations  must  be  done  using  this 
approach.  For  a  simple  search  through  T  =  40  points  with 
N  =  12,  5.5  billion  models  would  have  to  be  evaluated! 

This  approach  is  therefore  impractical  and  so  a 
sub-optimal  approach  which  yields  a  good  but  not 
necessarily  the  best  set  of  ' s  for  our  model  must  be 
used . 


If  we  view  this  problem  as  one  of  predicting  a 
dependent  variable,  Vi  f  rom  a  large  set  of  independent 
variables,  V^'s,  then  the  standard  linear  model  approaches 
may  be  applied.  In  a  statistical  sense,  independent 
variables  are  values  that  can  be  observed  but  not 
controlled  and  dependent  variables  are  affected  by  changes 
in  the  independent  variables.  Thus  the  value  of  dependent 
variables  is  said  to  depend  on  values  or  changes  in 
independent  variables.  The  linear  model  is  just  one 
approach  to  explaining  the  relationship  between 
independent  and  dependent  variables. 

In  most  linear  model  applications,  the  criterion  used 
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to  evaluate  the  fit  of  the  model  to  data  is  mean-square 

error.  It  is  desirable  in  most  cases  to  choose  a  model 

which  minimizes  the  mean-square  error.  The  problem  in  our 

case  is  to  choose  a  subset  of  \A  '  s  of  size  N  from  a  set  of 

1 s  of  size  T  such  that  the  linear  model  employing  those 

' s  produces  the  minimum  mean-square  error  when  compared 

to  all  other  possible  models  containing  N  V^'s.  This 

T 

cannot  be  done  without  examining  all  (N )  models  again 
however,  suboptimal  approaches  producing  a  very  good,  but 
not  necessarily  the  best  fit  are  available.  One  method 
employing  a  forward  selection  procedure  is  described  in 
the  following  sections. 

As  a  final  note  it  should  be  observed  that  we  are 
choosing  the  kernel  for  a  non-linear  N-gram-based  on  the 


fit  of  a  linear 

model 

to  sample  data. 

This  approach 

is 

admittedly  ad 

hoc 

and  is  chosen 

for  simplicity 

and 

computational  ease. 

Regardless,  it  is 

believed  that 

the 

kernel  chosen  by  this  approach  is  very  good  if  not  the 
best.  If  the  value  of  a  particular  is  important  to  the 
prediction  of  it  will  have  a  high  partial  correlation 
coefficient  when  it  is  examined  for  entry  into  the  model 
during  the  forward  selection  proceedure.  This  is  true  in 
all  cases  when  the  pixels  are  binary-valued  and  is  usually 
true  when  the  pixels  are  continuously  valued.  Still,  as 
we  are  not  considering  all  (£)  possible  sets  of  pixels  in 


the  neighborhood,  the  best  model  will  not  always  be  found 
[20]  . 


The  linear  model  which  we  use  to  determine  the  V^'s 
in  the  kernel  may  be  expressed  in  linear  regression  form 
as 

Yk  =  xj  Ek  k  =  1,2,  .  .  .  ,M 

where 


Yk  =  VN+l,k 


1 

Vl,k 

V2,k 


*-V, 


(3.4) 


N,k 


B  is  an  (N+l)xl  vector  of  unobservable  parameters  and 
is  an  unobservable  random  variable  such  that  E[ek]  =  0. 
The  sample  number  is  denoted  by  k  and  there  are  a  total  of 
M  samples.  We  can  also  define  matrices  X  and  Y  as 


X  = 


-  - 

— 

+T 

xi 

Y1 

X2 

Y2 

• 

Y  = 

• 

• 

T 

• 

.  V 

- 1 

S 

_ i 

(3.5) 


-*•  T  -1  T 

The  most  common  estimator  of  8  is  (X  X)  X  Y ,  the 


least-squares  estimator.  so  given  a  parent  texture  wnich 
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we  desire  to  simulate,  an  estimation  of  the  texture  model 
parameter  3,0,  may  be  obtained  yielding  the  generation 
model 

Y  =  xt  (3.6) 

It  is  very  important  to  note  that  the  estimation  of  the 
amount  of  noise  present  in  a  texture  is  usually  obtained 
by  measuring  the  amount  of  variation  in  the  parent  texture 
which  is  unexplained  by  our  model.  This  is  expressed 
numerically  in  the  mean-square  error.  Therefore,  as  in 
any  modelling  process,  any  shortcomings  or  inaccuracies  of 
the  model  will  appear  to  be  "noise"  (unexplained  variance) 
and  hence  the  amount  of  noise  will  increase. 

3.5  Correlation  and  Partial  Correlation 

We  may  define  the  mean  and  variance  of  a  variable  as 

Mv  =  EtV^  (3.7) 

and 

%  *  Et<v1-uVi>2]  <3,8) 

Similarly  we  may  define  the  covariance  of  two  variables 

and  V  as 
2 

Y„  „  =  E[ (V,-yv  ) (V,-yv  ) ] 

V1V2  1  1  ^  2 

(3.9) 


*  E<<vlv2>  -  “VjV,! 
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And  their  correlation  coefficient  as 


'V,  V 


12 


V1V2 


aV  aV 
V1  2 


(3.10) 


Using  the  above  definitions  we  may  then  define  the  partial 
correlation  coefficient  of  variables  V-^  and  V2  after  both 
have  been  adjusted  for  as 


w 


V. 


Pv1V2~(Pv1V3)(Pv2V3’ 

^  V1V32  ^'PV2V3 


(3.11) 


Each  of  the  above  parameters  has  a  corresponding  statistic 
(estimate  of  a  parameter  of  a  population  given  an 
observable  sample  of  that  population)  which  is  chosen  to 
meet  some  desirable  set  of  a  criteria  such  as  sufficiency, 
consistency  and  unbias  of  estimate  under  certain 
conditions.  Given  that  the  ^'s  are  samples  from  an  ith 
order  multivariate  normal  distribution  the  maximum 
likelihood  estimates  for  the  above  parameters  are 


(3.12) 


<vi,k-V2 


M 


(3.13) 
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E<Vl.k'Vl)(V2.k-V2> 


viv2 "  {[x;(vi,k-^>2][x;(v2,k-^»]2}i/2 


(3.14) 


r  —  ( r  )  ( r  ) 
V1V2  V1V3  V2V3 


'VlVV3 


Vl"rv  v 


(3.15) 


"  V1V3 


V2V3 


Second-order  partial  autocorrelation  may  be  found  in  a 
similar  manner  using 


^iVVXvVVV^ 


'VlVV3V4 


V1V4 ’ V3 


(3.16) 


V„V.  -v-3 
2  4  3 


Higher-order  partial  correlations  may  be  found  by 
extension  of  the  above  [21]. 


3.6  The  Forward  Selection  Procedure 


One  approach  to  finding  a  good  subset  of  independent 
variables  is  known  as  the  forward  selection  procedure. 
This  method  inserts  N  variables  into  the  model 
one-at-a-time.  The  order  of  insertion  is  determined  by 
using  the  partial-correlation  coefficient  as  a  measure  of 
the  importance  of  variables  not  yet  in  the  equation.  The 
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basic  p  idure  is  as  follows.  First  we  select  the 
most  correlated  with  VN+1  denoted  as  ,  and  we  find  the 

1  A 

first-order  linear  regression  equation  vN+i  =  aivi  +a2 • 
We  next  find  the  partial  correlation  coefficient  of  \T 
(for  all  j  /  i)  and  VN+^  (after  allowance  for  ). 


Mathematically  this  is  equivalent  to  finding  the 
correlation  between  the  residuals  from  the  regression 

A 

VN+^  =  a1Vi  +a 2  and  the  residuals  from  another  regression 
V.  =  b  V.  +b~  (which  we  have  not  actually  performed). 

]  i  i 

The  Vj  with  the  highest  partial  correlation  coefficient 

with  vN+i»v^  *  now  selected  and  a  second  equation 

V  =  c  V.  +c  V.  +c_  is  fitted  to  the  data.  This 
N+l  1  2  ±2  ^ 

process  continues.  After  V.  ,V.  ,...,V.  are  in  the 

X1  12  xq 

regression  the  partial  correlation  coefficients  are  the 
correlations  between  the  residuals  from  the  regression 

A 

VN+1  =  f (V^  ,V.  , . . . ,V^  )  and  the  residuals  from  a 

1  2  q 

regression  \T  =  f  .  (V  ^  )  ( j^q)  .  The  V.  with 

j  3  1  2  Xq  Xj 

the  highest  partial  correlation  coefficient  is  now 

selected  for  entry  into  the  linear  model.  The  process  is 

continued  until  N  ' s  are  entered  into  the  model. 


The  final  N  variables  chosen  by  a  forward  selection 


procedure  are  not  guaranteed  to  be  the  optimal  set  but 
given  the  logistics  of  the  selection  procedure,  the 
solution  obtained  is  usually  close  to  optimal. 


r 


One  o£  the  most  common  procedures  for  implementing 

the  forward  selection  process  numerically  utilizes 

Doolittle  decomposition  [22],  The  Doolittle  decomposition 

may  be  used  to  find  the  inverse  of  the  correlation  matrix, 

R  ,  and  the  estimates  of  linear  model  coefficients  as  each 
P 

variable  is  entered  in  the  model.  The  correlation  matrix 
merely  consists  of  the  set  of  correlation  coefficients 
rv  y  as  defined  by  Eq .  (3.14).  It  is  then  factored  into 

i  j 

the  product  of  two  triangular  matrices 


R 

P 


L  U 
P  P 


(3.17) 


where  L  is  lower  triangular  and  U  is  upper  triangular 
with  ones  on  the  diagonal.  Partial  autocorrelation 
coefficients  may  be  obtained  easily  from  elements  of  this 
matrix  during  its  decomposition  at  each  step.  For  further 
details  and  examples  on  the  forward  selection  procedure  of 
the  decomposition  process  see  Beyer  [18]  and  Draper  and 
Smith  [20] . 


3.7  Statistics  Collection 


In  practice,  N,  the  number  of  pixels  in  the 
generation  kernel,  is  often  chosen  by  computational  limits 
imposed  by  finite  processing  capability  or  finite  computer 
memory.  The  idea  of  parsimony  would  also  urge  the 
selection  of  the  smallest  N  possible.  In  the  most  general 
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texture  synthesis  (stochastic)  model  which  utilizes 

N 

N-grams  we  find  that  as  many  as  g  storage  locations  are 
required  in  order  to  collect  the  data  needed  to  synthesize 
a  texture  from  a  parent  where  g  is  the  number  of  grey 
levels  in  the  texture.  This  approach  calls  for  N  to  be 
less  than  17  for  g  =  2  (a  binary  image)  if  we  have  only 
2^  =  131,072  storage  locations  in  memory.  Approaches  to 
"stretch"  this  limitation  have  been  investigated  and  will 
be  discussed  in  a  latter  section.  If  we  have  8  grey 
levels  then  N  must  be  less  than  or  equal  to  5  as 

5 

8  =  32,768.  Thus,  in  synthesizing  textures  using  an 
N-gram  approach,  processor  storage  capability  is  the  major 
limiting  factor. 

Determining  which  pixels  will  be  included  in  the 
generation  kernel  requires  an  estimate  of  the  linear  model 
defined  in  Eq .  (3.4).  To  do  this  the  X  and  Y  matrices  of 
Eq .  (3.5)  may  be  used  or  the  correlation  matrix  of  the 
kernel  points  must  be  estimated  using  a  parent  texture. 


The 

elements 

o  f 

the 

*k 

vector  of 

Eq . 

(3.4)  are 

obta i ned 

by  passing 

the 

kernel  window  over 

the  sample 

parent  texture  and 

record ing 

the 

pixel  value 

correspond ing 

to  the 

po s  i  t  ion 

o  f 

each  V  ^ 

i  n 

the  window. 

The 

kernel  is 

assumed  to  be  completely  within  the  boundary  of  the  parent 
texture.  A  sample  may  be  taken  at  all  possible  positions 


in  the  parent  texture  or  a  random  sample  of  points  may  be 
chosen  if  the  parent  texture  is  very  large  to  reduce  the 
number  of  computations  required.  In  actual  practice,  no  X 
matrix  is  ever  formed.  To  obtain  the  correlation  matrix 
of  the  sample  set,  R,  (and  from  that  6)  we  need  only  the 
surn  of  squares  cross  products,  and  V  ^ ,  over  the  sample 
set  according  to  Eq .  (3.14). 

The  elements  of  the  correlation  matrix  for  many 
kernel  patterns  often  contain  redundancies.  For  example 
the  spatial  relationship  between  the  pair  V-,  and  V2  and 
the  pair  V3  and  V5  in  Fig.  3.4  is  the  same.  Estimating 
the  correlation  for  these  two  pairs  of  points  from  a 
sample  will  yield  nearly  equal  values  if  the  sample  size 
is  large  and  the  overlap  of  samples  used  to  estimate  the 
correlation  values  is  large.  In  short, 

rV.V„  =  rV-.V.  (3.18) 
12  3  4 

In  more  mathematical  terms,  let  Ifn^n^  denote  the 
random  texture  field  where  n  ^  and  n^  are  integers 
representing  the  coordinates  of  points  in  our  sampled  and 
quantized  image.  Let  n  be  the  vector  of  coordinates 
(n^,n^).  From  Section  2.1,  second-order  or  2-gram 
statistics  are  given  by  the  set  of  joint  density  functions 

P+  -*(vi'v-j)  (3.19) 
n  ,m  J 
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for  all  possible  vectors  n  and  in,  where  and  V2  are  the 
pixel  values  of  the  random  variables  I  (n)  and  I  (in)  , 
respectively.  If  the  random  field  is  homogeneous,  that 
is,  invariant  through  translations  then 


P  =  P 

-+  -y  -v 

n,m  n-m,  0 


from  Eq .  (2.2). 


(3.20) 


If  we  assume  that  our  two-dimensional  texture  is 
homogeneous  and  stationary  we  might  propose  that 
Eq .  (3.20)  is  true  by  definition.  A  word  of  caution  is 

necessary  at  this  point.  Estimating  the  elements  of  the 
covariance  or  correlation  matrix  using  the  assumption  of 
stationarity  and  homogeneity  can  yield  correlation 
matrices  having  negative  determinants,  in  violation  of  the 
fact  that  non-singular  covariance  matrices  must  be 
symmetric  positive  definite.  The  violation  occurs  because 
the  sample  is  not  homogeneous  and  stationary.  Therefore 
the  type  of  assumption  expressed  in  Eq .  (3.20)  should  only 
be  used  for  estimation  when  sample  sizes  are  very  large 
and  are  known  to  be  homogeneous  and  stationary. 

The  assumption  of  Eq .  (3.20)  can  be  very  powerful  in 
a  computational  sense  for  large  samples  as  the  calculation 

A 

of  covariance  (or  correlation)  matrices  can  be  time 
consuming.  The  complete  covariance  matrix  Cor  a 
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contiguous  10  x  10  (100  element)  kernel  contains  10,000 


entries.  Of  these,  4050  are  redundant  by  simple  symmetry 


as  /  jM  V  =  /  jV  V  which  implies 
l  g  j  l 


r  =  r 
V  •  V  ■  V  •  V  • 

i  j  v  g  i 


Computing  all  cross  products  and  sums  to  estimate  this 
matrix  from  a  512  x  512  image  requires  over  1.025  billion 
multiplications  and  1.05  billion  additions.  Utilizing  the 
concept  of  homogeneity,  this  may  be  reduced  to  0.05 
billion  mul t i pi icat ions  and  0.075  billion  additions  which 
is  significantly  lower. 


Once  the  covariance  matrix  of  a  kernel  is  determined, 
the  linear  model  containing  the  V^'s  of  the  kernel  may  be 
obtained.  It  is  advisable  to  make  the  kernel  large  as 
more  texture  information  is  contained  over  large 
distances.  In  many  cases,  there  are  relationships  between 
pixels  separated  by  great  distances  especially  if  the 
texture  is  coarse  or  highly  regular  (periodic).  In 
practice,  however,  large  matrices  may  be  numerically 
ill-conditioned  during  a  decomposition  or  inversion 
process  and  so  they  must  be  avoided.  Also  computer 
storage  limitations  must  be  considered.  In  this  study,  no 
matrix  larger  than  100x100  was  decomposed  or  inverted  for 
these  reasons.  This  constraint  required  a  mul t i pi  e-pass 
approach . 


First,  100  V^'s  (usually  closest  to  the  kernel  eye. 
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V  ,)  were  chosen  for  examination.  They  were  entered  into 

the  linear  model  in  a  forward  selection  manner  until  it 

was  determined  that  entry  of  additional  V^'s  would  be 

insignificant.  Insignificance  is  indicated  both 

statistically  by  considering  the  reduction  in  the  sum  of 

squares  due  to  the  entry  of  a  variable  into  the  model  [20] 

and  also  by  the  value  of  6  which  approaches  zero  so  the 

i 

variable  V.  becomes  insignificant.  On  the  second  pass, 
variables  not  tried  in  the  first  pass  are  examined  and 
tested  for  possible  entry  into  the  model.  Again,  those 
which  are  significant  are  kept  and  the  others  are 
discarded.  On  the  next  pass,  any  variables  not  examined 
in  the  first  two  passes  may  be  examined.  The  process 
continues  until  all  V^'s  have  been  tested  for  possible 
entry  into  the  linear  model  at  least  once. 

This  multiple-pass  process,  as  the  forward  selection 
procedure  which  it  employs,  is  not  guaranteed  to  produce 
the  best  model  but  will  provide  an  excellent  model 
nevertheless.  The  final  model  will  also  be  parsimonious. 
This  will  reduce  the  number  of  parameters  in  our  model  and 
aid  in  the  efficiency  of  the  generation  process. 

3.8  Results 

Once  the  points  for  the  kernel  are  chosen  based  on 
the  linear  model  derived  using  the  methods  described  in 
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the  previous  sections,  estimates  o£  the  generation 

parameters  tor  the  texture  are  obtained  using  concepts 

discussed  in  section  3.2.  Practical  considerations 

require  us  to  limit  N ,  the  number  of  pixels  in  the  kernel, 

to  12  to  13  depending  on  the  processor  storage  available 
N 

as  2  values  must  be  stored.  These  G's  are  then  used  to 
generate  each  pixel  along  a  row,  row  by  row  until  a 
complete  two-dimensional  texture  is  obtained.  For  eacn 
pixel  the  appropriate  generation  parameter  estimate  is 
found  and  a  uniformly-distributed  pse udo- random  variable 
is  generated.  Based  on  these  two  values,  a  black  pixel 
(0)  or  white  pixel  (1)  is  generated. 


In  practice,  not  all  of  the  generation  parameters  may 
be  estimated  when  N  is  large  because  all  possible  patterns 

of  vi»  V2  '  *  *  *  '  vn'  VN+1  ma^  not  *3e  Present  if  the  sample 
image  or  there  may  be  few  of  them.  Smaller  samples  can 
cause  inaccurate  estimation  of  the  G's  as  the  variance  of 
our  estimate  is  larger  and  therefore  the  expected  error  of 
our  estimate  is  larger  than  would  be  expected  with  a 
larger  sample  size.  In  these  cases  it  is  important  to  sum 
over  the  least  significant  kernel  elements  and  estimate 

by  G„  (V.,V.,,,...,V„).  In  our 


°V  <VV2 . V 

N+l 

study,  this  was  done  if  the  sample  size  to  compute 


-vN+1(vi’Vi . V- 


G  (V  , V  , . . . , V  )  was  less  than  10.  The  variable  i  is 
VN+1  1  2  N 

increased  until  this  condition  is  met. 
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'i-..  uure  simulations  using  this  method  are  shown  in 
Figs.  3. 5 (b) -3 . 14 (b) .  Visually,  the  results  are  very 
good.  As  the  estimated  texture  generation  parameters  are 
approximated  using  statistics  gathered  from  the  full 
parent  texture,  non-homogeneity  in  the  parent  texture  will 
cause  an  "average"  texture  to  be  synthesized.  The 
simulation  of  straw  (Fig.  3.7)  is  poor  because  of  its 
non-homogeneous  nature  (specifically  the  directionality  of 
the  stalks  in  different  parts  of  the  image)  and  exhibition 
of  detail  (specifically  individual  non-conforming  single 
stalks).  A  similar  observation  may  be  made  with  respect 
to  the  parent  textures  of  grass  and  water  but  in  these 
textures  the  non-homogeneity  is  not  so  pronounced.  As  we 
are  attempting  to  synthesize  textures  and  not  merely 
"image  code"  the  parent  textures,  details  and 
non-homogeneities  will  be  lost  in  the  synthesis  process. 

The  bark  texture  is  among  the  most  difficult  to 
simulate  due  to  its  very  unusual  macro-structure.  Still, 
the  N-gram  simulation  looks  remarkably  similar  to  the 
original  when  windowed  regions  20  to  40  pixels  square  are 
obse  rved . 

* 

The  kernels  used  to  generate  these  N-gram  simulations 
are  shown  in  Table  3.1.  A  clear-cut  relationship  between 
the  textures  and  the  kernels  found  using  the  linear  model 
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Figure  3.5  Grass 


(a)  Original  Texture 


(b)  N-gram  Simulation 


(c)  Linear  Model  Simulation 
Figure  3.9  Leather 
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Original  Texture 


(b)  N-gram  Simulation 


(c)  Linear  Model  Simulation 


Figure  3.12  Wood 


(a)  Original  Texture  (b)  N-gram  Simulation 


(c)  Linear  Model  Simulation 


Figure  3.13  Pigskin 


Table  3.1  TWO-DIMENSIONAL  N-GRAM  TEXTURE 
SIMULATION  KERNELS 


GRASS  BARK 
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approach  is  not  always  evident.  Some  kernels,  such  as 

those  for  wool  and  water,  seem  to  have  a  vertical  or 

horizontal  structure  possibly  resulting  from  the  structure 
of  their  corresponding  textures  while  others,  such  as  the 
kernels  for  wood  and  raffia,  do  not.  Still,  the 

aon-cont iguous  kernels  produce  excellent  results. 

3.9  Extension  of  N-Gr am  Model 

The  N-gram  model  may  be  extended  beyond  the  fixed  w 

set  by  limited  processor  storage.  The  requirement  that  G 

must  be  estimated  based  on  no  fewer  than  q  (which  was  set 

to  10  in  our  study)  samples  already  reduces  the  storage  to 

a  maximum  of  M/q  non- red undant  parameters  where  M  is  the 

total  number  of  samples  in  the  parent  texture  image.  For 

5 

a  512x512  image,  M  is  approximately  2.5x10  .  Given 

q  =  10,  this  implies  that  a  maximum  of  25,000  generation 
parameters  estimates  must  be  stored.  3y  increasing  q, 
further  reduction  is  possible.  Meanwhile  N  is  permitted 
to  increase  without  bound  to  sample  size  limitations.  For 
frequently  occurring  patterns,  larger  N's  allow  increasing 
dependency  over  larger  distances.  This  is  most  desirable 
in  regular  and  coarse-structured  textures. 

Searching  for  the  proper  generation  parameter  at  each 
step  in  this  type  of  process  is  complex.  In  the  earlier 
N-gram  storage,  each  pattern  of  0's  and  l's  actually  forms 


x 
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a  binary  address  to  the  location  in  memory  of  the  desired 
G.  In  the  extended  case  a  sort,  search,  or  a  series  of 
comparisons  along  with  some  intelligent  preprocessing  is 
required.  Efficiency  is  reduced. 

To  illustrate  the  effect  of  model  extension  the 
texture  raffia  was  used.  Figure  3.14(a)  shows  the 
original  and  Fig.  3.14(b)  is  the  synthesis  obtained  using 
a  kernel  containing  N  =  14  pixels.  Figure  3.14(c)  is  the 
linear  model  synthesis.  Figures  3. 15  (a)  , ( b) , ( c)  were 
obtained  using  three  different  texture  kernels  with  N  =  22 


points . 

Far 

more  structure 

in 

these  extended 

model 

versions 

is 

apparent.  This 

is 

expected  as  at 

each 

generation  step,  the  next  pixel  is  allowed  to  depend  on 
pixels  further  from  it.  As  the  pixels  in  the  kernel 
become  more  widely  spaced  the  synthesis  becomes  more 
structured  but  small,  local  regions  often  become  more 
distorted  and  less  r a f f i a- looki ng  because  the  information 
used  in  the  synthesis  process  is  more  global  than  local. 

3.10  Linear  Model  Generation  of  Binary  Textures 

The  process  of  choosing  the  V  '  s  to  be  present  in  the 

l 

texture  generation  kernel  described  in  earlier  sections  of 
this  chapter  actually  yields  a  simple  linear  model  which 
can  also  be  used  to  generate  binary  textures.  The  model 
which  results  from  the  determination  of  the  generation 
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kernel  may  be  expressed  in  equation  form  as 


VN+l,k  "  8lVl,k+  82V2,k+"-  +  eNVN,k+  V  £k 


(3.21) 


or  more  simply  as 


VN+1  '  BlV  82V2+---+  W  V  e  ’ 


(3.22) 


Once  the  estimates  of  the  B^'s  are  known,  a  pixel  VN+1  may 
be  calculated  from  a  set  of  given  values  plus  an  error 
e.  In  one-dimensional  analysis  this  is  sometimes  known  as 
the  autoregressive  time  series  model  [24J.  For  binary 
a  value  of  VN+^  will  be  produced  which  is  non-binary.  To 
generate  binary  data  using  this  model  will  therefore 
require  quantization. 


In  the  N-gram  approach  to  texture  simulation,  the 
randomness  of  the  texture  is  induced  by  the  generation  of 
a  uniformly-distributed  pseudo- random  variable  during  the 
generation  process.  The  comparison  of  this  value  with  the 

A 

estimate  of  the  generation  parameter,  G,  yields  the  next 
binary  pixel.  A  similar  type  of  randomness  must  occur  in 
the  generation  of  binary  textures  using  the  linear  model 
of  Eq.  (3.22).  This  randomness  is  expressed  in  the  model 
in  the  error  term  e. 


We  can  obtain  an  estimate  of  the  distribution  of  e  in 
the  same  manner  as  we  estimate  the  0's  of  the  model.  This 


05 


may  be  done  by  applying  the  model  to  the  sample  data  from 
which  it  was  derived  and  observing  the  errors.  That  is, 
the  linear  aiodel  kernel  is  passed  over  the  parent  texture 

A 

image  and  at  each  point  a  V  is  calculated  based  on 

N+l 

A 

Eq.  (3.22)  without  the  error  term.  Then  V  -V  is 

N+l  N+l 

calculated  where  V  is  the  actual  value  of  V  in  the 

N+l  N+l 

parent  texture.  The  histogram  of  the  values  can  be  used 
to  estimate  the  distribution  of  £.  As  one  step  further  we 
could  assume  that  has  some  known  distribution  such  as 
Gaussian  or  normal,  and  merely  estimate  the  parameters 
necessary  to  define  this  distribution.  In  the  normal 
distribution  case,  only  the  standard  deviation  (or 

variance)  of  needs  to  be  estimated.  The  mean  of  e  is 

zero  in  the  linear  model,  least-squares  distribution. 

Our  generation  process  then  consists  of  the 

calculation  of  ^  ^  V.;  +  8Q  to  which  we  add  a  random, 

normally-distributed  error  term  e  and  this  value  is  then 

quantized  to  0  or  1  based  on  comparison  with  0.5.  Results 

using  this  generation  method  are  shown  in  Figs.  3.5(c) 

through  3.14(c).  In  these  figures,  N  was  allowed  to  be  as 

N 

large  as  70  as  only  N  coefficients  (not  2  )  need  to  be 

stored  along  with  ,  the  estimate  of  the  error  standard 

deviation. 

The  kernels  used  in  the  linear  model  simulations  are 
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p 


shown  in  Table  3.2.  The  linear  model  simulations  are 
slightly  inferior  to  the  N-gram  simulations  but  the 
degradation  is  far  less  than  we  would  expect  from  such  a 
massive  compression  of  information  (which  is  approx imately 
2  to  70).  The  results  were  good  enough  to  encourage  the 
application  of  the  linear  model  to  continuous-tone 
textures . 
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Table  3.2  TWO-DIMENSIONAL  BINARY  LINEAR  MODEL 
TEXTURE  GENERATION  KERNELS 


GRASS 


BARK 

cP  a 
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□  □ 


CHAPTER  4 


ALGEBRAIC  RECONSTRUCTION  TECHNIQUE  MODEL 

4. 1  Introduction 

Julesz's  conjecture  [8]  that  second-order  statistics 
are  sufficient  for  human  visual  texture  d isc r im ina t ion 
provides  a  useful  estimate  of  the  amount  of  information 
necessary  to  reconstruct  a  texture  field.  Although 
counterexamples  to  that  conjecture  have  recently  been 
found  [4,25],  and  are  shown  in  Chapter  2,  it  is  a  good 
first-order  approximation.  Examples  of  the  use  of  that 
upper  bound  for  texture  analysis  can  be  found  in  [14,26]. 
Therefore  it  is  very  tempting  to  use  it  for  synthesizing 
natural  texture  fields.  That  is,  we  may  attempt  to 
simulate  textures  based  on  generation  parameters  estimated 
from  2-gram  statistics  over  various  distances.  This 
approach  requires  fewer  statistics  to  be  collected  from  a 
parent  texture  as  only  2-grams  versus  N-grams  are 
collected . 

This  chapter  illustrates  that  we  must  "invent" 
higher-order  statistics  to  use  the  Markov  generation 
coefficient  approach  for  texture  synthesis  if  we  limit  our 
knowledge  of  the  original  random  field  to  second-order 
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statis  ts.  We  will  demonstrate  how  this  can  be  done 
using  Algebraic  Reconstruction  Techniques  (ART). 

4.2  Problem  Definition 

The  stochastic  approach  toward  texture  analysis 
considers  texture  fields  as  samples  of  2-D  stochastic 
fields.  Second-order  statistics  are  given  by  the  set  of 
second-order  joint  density  functions  (2-grams),  P(V^,Vj). 
As  in  Chapter  3,  we  also  assume  for  simplicity  that  the 
random  field  is  homogeneous,  that  is,  invariant  through 
translations.  Given  the  set  of  joint  density  functions, 
(P(VifVj),  for  all  spatial  relationships  between  and 
V j ,  our  problem  is  to  synthesize  a  texture  field  with  the 
same  second-order  statistics. 

There  are  many  ways  of  carrying  out  such  a  synthesis. 


We  use 

a 

television 

raster 

type  of  scanning  (left 

to 

right , 

top 

to  bottom) 

and 

the  generation  kernel 

of 

Fig.  3. 

1(b)  . 

Even  if  we 

assume 

finite  memory,  namely  that 

intensity  at  point  VN+j  depends  only  on  intensities  at 
points  located  in  some  finite  neighborhood,  we  see  that 
second-order  statistics  are  not  sufficient  to  generate  the 
process.  Indeed,  assuming  a  memory  of  order  N,  the 
intensity  at  V  is  computed  using  the  conditional 
probability  or  generation  coefficient 

Gy  (V^jVji  •  •  •  iVjj)  =  ‘  ’ '  (4.1) 
N+ 1 
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which  involves  (N+l ) th-order  joint  density  functions. 
Therefore  if  N  is  larger  than  1  and  we  are  given  only 


second-order  statistics,  P(\T,\T),  we  have  to  "invent" 
higher-order  densities  PfV^V^,  ...,V  ^)  .  Mathematically, 

the  problem  can  then  be  stated  as  follows:  given  N  random 
variables  vi'*,-'Vn  SLJch  that  for  every  pair  (vi'vN)» 
l^i,j<N  and  i  ^  j ,  we  know  the  joint  density  function 
P(Vi,Vj),  find  a  function  P(V^,...,VN)  which  satisfies 


P (V1 , V2 » • • • / VN)  >  0  for  all  V1'---'VN  (4.2) 


£••£  £•••£  £••£-, . V- 


V.  V.  ,  V.  ,  ,  V  .  .  V.  ,  V.T 

1  1-1  1+1  ]-l  3+1  N 


* Vj  f • • •  < - 
(4.3) 


P(Vi,Vj)  . 

Here  the  P (V  ,V  )'s  are  sometimes  called  the  marginals  of 
i  ] 

P (V  ,...,V  ).  Assuming  quantization  with  g  levels,  the  gN 
IN 

unknowns  P  (V  ,  .  .  .  ,V  )  can  be  stacked  as  a  vector  and 

I  N 

conditions  (4.2)  and  (4.3)  correspond  to  a  linear 
programming  problem: 


IAp  =  rn  (4.4) 

p  >  0  (4.5) 

where  vector  is  obtained  from  the  functions  P(V^,V^)  and 
matrix  A  of  size  ((^)9^)X9N  contains  only  ones  and  zeroes. 
For  any  reasonable  values  of  N  and  g,  this  is  a  set  of 

linear  equations  and  inequalities  of  fairly  large 
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dimension  and  the  usual  solution  techniques  such  as  the 
simplex  method  [4]  become  very  limited. 

4.3  Solution  Through  Algebraic  Reconstruction 

Techniques  (ART) 

The  ART  algorithm  was  introduced  by  Gordon,  Bender 
and  Herman  [27]  for  solving  the  problem  of 
three-dimensional  reconstruction  from  .projections  in 
electron  microscopy  and  radiology.  This  is  a 
deconvolution  problem  in  which  a  function  in  a 
higher-dimensional  space  is  estimated  from  its 
experimentally  measured  projections  in  a  lower-dimensional 
space.  For  an  excellent  review  of  those  techniques  see 
Gordon  [ 28 ] . 

The  problem  stated  in  Eqs.  (4.2)  and  (4.3)  or  (4.4) 
and  (4.5)  is  precisely  of  this  form,  where  the  projections 
are  the  second-order  joint  density  functions.  ART  is 
therefore  directly  applicable.  The  simple  intuitive 
interpretation  is  that  each  projected  density  is  thrown 
back  across  the  higher-dimensional  region  from  whence  it 
came,  with  repeated  corrections  to  bring  each  projection 
of  the  estimate  into  agreement  with  the  corresponding 
measured  projection. 

Formally,  we  use  an  iterative  scheme  defined  by 
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P(q+1)  (V,  , . . .  ,v  >  =  P(g)  (V,  ,  .  .  .  ,v„)+tc  (q)  (V, , . . .  ,V) 


(4.6) 


For  all  values  of  V^,...,VN,  for  q  =  0,1, 


where  the  correction  term  c^)  (Vj_  #  . . .  ,Vjj  )  is  given  by 


c  ^  (V1 , • • • »VN)  - 
N-l  N 


(4.7) 


(2J  i=l  j=i+l 


and  P(Vi,Vj)  is  the  actual  marginal  measured,  for  example, 

%  (q) 

from  an  original  texture  field.  Here  P  (V  ,V  )  is  the 

i  j 

marginal  corresponding  to  the  reconstructed  density  at 
iteration  q. 


We  may  express  this  in  words  as  follows.  The 

iterative  process  may  be  started  with  all  reconstruction 

,  ^ (0)  ,  .  n 

elements  set  to  a  constant  (P  (V  ,...,v  )  =l/g  for  all 

1  N 

(V  ,...,V  ).  in  each  iteration  the  sum  of  the  differences 
1  N 

between  the  actual  and  the  reconstructed  marginals  is 

N— 2 

computed  and  evenly  divided  amongst  the  g 
reconstruction  elements.  If  the  correction  is  negative, 
it  may  happen  that  the  calculated  density  becomes  negative 
at  some  points.  This  problem  can  be  alleviated  by  using  a 
modified  iteration  scheme  defined  by 


p(q+1>  (Vr  .  .  .  ,VN)  =  MAxjo,P<q)  (Vj^  ..  .  ,VN)+tc<q)  (Vj,..  ..VN)} 

Mq+1)  (4-8> 

therefore  guaranteeing  P  0  (constrained  ART  [28]). 


It  is  of  course  necessary  to  determine  when  an  iterative 
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algorithm  has  converged  to  a  solution  which  is  optimal 
according  to  some  criterion.  This  is  in  turn  related  to 


the  problem  of  finding  the  optimal  value  t  of  t.  Various 
criteria  for  convergence  have  been  devised  [28],  For 
simplicity,  we  chose  the  mean-square  error 


£  =  ||<£ 
q 


(q) 


2 

2 


i  -»■  'v  (q)  |  ,  2 

I"1"®  I  I 2 


(4.9) 


between  the  measured  and  calculated  marginals  where  | | • |  L 
is  the  usual  euclidean  no/m  and  in  is  defined  in  Eq.  (4.4). 


To  derive  the  optimal  step  size,  t,  for  each 
iteration  we  rewrite  Eq .  (4.6)  in  vector  form  as 


=  £<q>+ 


(4.10) 


Multiplying  both  sides  of  Eqs.  (4.10)  with  matrix  A 
(Eq.  (4.4)  we  obtain 


£(q+i)  =  ^ (q) 


m  +  td 


(q) 


(4.11) 


where 


2<q>  „  A  3<q) 


(4.12) 


and  subtracting  the  actual  marginal  vector  m  from  both 
sides  of  Eq .  (4.11)  yields 

||S(q+l,||2  „  ||j(q).  ta<q)||2 

=  t2 1 id <CI) m |-2ta -e (C1)  + 1 1 4 ''J’ 1 1 1  !4‘“) 


Therefore  the  error  e  ,  at  iteration  o+l  is  minimized  for 

q  +  1 
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A 


(4.i4: 


=  gtq>.a(q> 
tq  "  ll3(q,M* 

where  •  denotes  the  inner  product. 

A  dual  approach  also  explored  is  based  upon  the 

analysis  of  Eq .  (4.3)  in  the  Fourier  domain.  It  car 
easily  be  shown  that  the  initial  problem  stated  by 
Eqs.  (4.2)  and  (4.3)  is  equivalent  to  an  interpolation 

problem  in  the  Fourier  domain.  The  major  drawback  of  this 

approach  is  the  difficulty  to  ensure  the  positivity  of  the 
inverse  Fourier  transform  of  the  interpolated  function. 
Therefore  this  method  was  not  pursued  even  though  it  may 
be  the  case  that  "good"  interpolating  functions  will 
alleviate  that  problem. 

The  basic  philosophy  of  the  two  approaches  just 

discussed  is  that  Nth-order  joint  density  functions  are 

"invented"  to  satisfy  exactly  the  constraints  stated  in 

Eq.  (4.3).  Their  obvious  disadvantage  is  the  high 

dimensionality  (gN  for  an  Nth-order  joint  density 

function)  of  the  data  that  is  to  be  stored  compared  with 

N  2 

the  usually  lower  dimensionality  ( (^ ) g  for  the 
second-order  joint  density  functions)  of  the  data  that  is 
effectively  used. 
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4.4  Results  and  Conclusions 


The  iterative  process  defined  in  Eq.  (4.6)  may  be 
halted  at  any  number  of  iterations,  q,  and  a  texture  may 

-4- 

be  generated  using  the  value  of  p  at  that  point.  However, 
it  should  be  kept  in  mind  that  the  success  of  a  texture 
synthesis  depends  on  making  the  error  c ^  as  small  as 
possible  and  that  the  texture  generation  process  is 

sensitive  to  this  error.  It  has  also  been  found  by 

% 

experimentation  that  the  p  contains  many  values  which  are 
set  to  zero  by  implementation  of  constrained  ART.  This 
tends  to  cause  the  Markov  texture  generating  process  to 
become  absorbing,  which  causes  patches  of  white  and  black 
or  streaks  and  lines  to  be  generated.  This  is  eliminated 
by  setting  those  values  which  are  zero  to  a  small  non-zero 
value,  6  ,  in  the  generation  process. 

Using  the  above  concepts,  texture  simulations  of  the 
binary  textures  water  (Fig.  4.1(a))  and  raffia 
(Fig.  4.2(a))  were  generated  (Figs.  4. 1  (b)  ,4. 2  (b)  )  . 
Textures  similar  to  those  in  Chapter  3  employing  actual 
Nth-order  statistics  (N=14)  were  also  generated 
(Figs.  4 . 1  (c)  , 4 . 2  (c) )  and  are  included  here  for  reference. 
The  in  the  generation  kernels  for  the  ART  model  and  the 
N-gram  model  were  chosen  using  the  ideas  described  earlier 
in  Chapter  3. 
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(b)  ART-generated  Simulation 


At  the  onset  of  investigation  of  the  Algebraic 
Reconstruction  Model  it  was  hoped  that  this  approach  would 
be  useful  for  generating  continuous-tone  textures. 
However,  using  it  to  generate  binary  textures  revealed 
that  the  convergence  for  the  iterative  process  was  very 
slow  even  with  the  optimal  step  size  t.  Each  iteration 
required  much  computation  and  the  storage  required  for  the 
N-grams  was  large.  More  than  eight  hours  of  CPU  time  on  a 
DEC  KL-10  were  required  to  execute  the  large  number  of 
iterations  required  for  a  visually  pleasing  solution. 
Generating  textures  based  on  estimated  N-grams  which  is 
detailed  in  Chapter  3  is  probably  more  efficient  and  less 
complex  computationally.  This  work  did  lead  to  other 
texture  simulation  models  (see  Faugeras  [29]). 

The  results  using  algebraic  reconstruction  are  nearly 
equivalent  to  the  results  using  complete  N-grams  as 
second-order  statistics  collected  from  binary  textures 
contain  a  great  deal  of  information.  Close  examination  of 
the  textures  generated  using  algebraic  reconstruction  on  a 
high-resolution  display  device  reveals  a  high  random  noise 
level.  The  computational  requirements  and  the  final  noise 
level  indicate  that  the  N-gram  method  of  generating 
textures  is  much  less  complex  and  yields  better  visual 
results  than  the  method  employing  algebraic 
reconstruction. 
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CHAPTER  5 


CONTINUOUS-TONE  LINEAR  TEXTURE  MODEL 

5.1  Introduction 

In  this  chapter,  the  concept  of  using  a  linear  model 
for  generating  binary  textures  that  was  briefly  discussed 
in  Chapter  3  is  extended  to  multi-grey-level  (256-level) 
or  continuous- tone  textures.  Pictorial  results  of  the 
application  of  this  model  to  simulation  of  a  variety  of 
textures  are  presented. 

The  application  of  the  linear  or  autoregressive  model 
to  time  series  processes  has  been  extensive.  These 
applications,  which  range  from  weather  forecasting  to 
stock  market  predictions,  primarily  utilize  the  models  and 
concepts  introduced  by  Box  and  Jenkins  [24],  Many  authors 
have  expanded  and  elaborated  on  their  approaches 
[30-38,62].  Some  researchers  have  applied  these  ideas  to 
texture  simulation,  in  fact,  the  Box-Jenkins 
autoregressive  model  is  one  of  the  very  few  approaches  to 
texture  synthesis  presented  thus  far. 

McCormick  and  Jaya ramamur thy  [2]  were  perhaps  the 
first  to  make  a  notable  attempt  to  simulate  natural 
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textures  using  this  approach.  Their  work  consisted  of  a 
discussion  of  the  Box  and  Jenkins  autoregressive  (AR) , 
moving  average  (MA)  and  autoregressive  integrated  moving 
average  (ARIMA)  models  including  estimation  of  model 
parameters  and  adequacy  of  model  fit.  (These  terms  are 
later  defined  in  this  chapter.)  A  very  simple  model  was 
then  used  to  simulate  two  very  similar  textures  which 
closely  resemble  the  wood  texture  of  this  study  by  filling 
in  the  holes  of  a  parent  texture  using  the  derived  model. 
Only  two  textures,  both  exhibiting  a  wood-grain-like 
structure,  were  used.  Similar  work  was  done  later  by  Tou, 
Kao  and  Chang  [11].  Unfortunately,  the  results  of  their 
simulation  of  these  textures  were  displayed  using  a 
printout  of  Chinese  cYiaracters  and  so  the  degree  of 
success  of  their  method  is  unclear.  The  appearance  of 
texture  synthesis  results  on  a  computer  printout  will 
confuse  most  observers  unaccustomed  to  such  crude  image 
displays.  The  models  were  again  very  simple  and  contained 
no  more  than  three  terms  in  the  linear  model  summation. 
Deguchi  and  Morishita  [12]  attempted  to  use  the  linear 
model  to  segment  and  partition  textures.  Their  approach 
was  only  partially  successful. 

In  the  above  simulation  attempts,  the  models  used 
were  simple.  The  process  of  collecting  statistics  and 
estimating  parameters  is  complex.  In  some  cases,  previous 
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authors  attempted  to  use  the  complex  Box  and  Jenkins  ARIMA 
model  which  leads  to  difficult  model  parameter  estimation 
if  the  number  of  model  elements  is  greater  than  two  or 
three . 


In  our  study,  the  simpler  autoregressive  model  is 
used  and  is  allowed  to  contain  a  large  number  of 
parameters.  This  is  possible  using  the  assumption  of 
homogeneity  (stationarity)  combined  with  the  forward 
selection  process  of  choosing  non-contiguous  generation 
kernels  as  described  in  Chapter  3.  These  models  are 
extended  further  by  allowing  second-order  autoregressive 
models  and  non-stationary  noise.  Results  of  texture 
simulations  using  these  models  are  included  in  this 
chapter . 

5.2  The  Linear  Autoregressive  Model 

In  Chapter  3  the  linear  autoregressive  model,  used  to 
determine  the  elements  of  the  generating  kernel,  was 
expressed  as 

Yk  =  *k^  +  ek  k  =  1, — ,M  (5.1) 

where 

Y  =  V 
k  N+l ,k 
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and 


Xk  = 


1 1  k 
2  ,  k 


L  v 


N,k 


Here  8  is  an  (N+l)xl  vector  of  unobservable  parameters  and 

e  is  an  unobservable  random  variable  such  that  E[e, }  =  0. 
k  k 

The  sample  number  (index)  is  denoted  by  k  and  M  is  the 
total  number  of  observations.  We  can  also  define  the 
vectors  Y  and  t  and  the  matrix  X  by 


r  ±T 

r 

r  _ 

X1 

Y1 

ei 

X" 

Y„ 

£ 

2 

2 

2 

-> 

• 

Y  = 

* 

e  = 

• 

• 

• 

• 

L  xM 

•  ym- 

•eM" 

and  our  model  may  be  expressed  as 

Y  =  X$  +  e  (5.3) 

In  equation  form,  dropping  the  k  subscript,  the  model 
becomes 
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Sums  and  sums  of  squares  leading  to  the  calculation  of  the 
correlation  or  covariance  matrix  of  the  parent  texture  are 
obtained  by  passing  any  chosen  generation  kernel  pattern 
over  the  texture.  From  this  matrix,  the  least-squares 
parameter  estimate  of  6  is  obtained.  The  multiple-pass 
forward  selection  process  described  in  Chapter  3  leads  to 
a  final  linear  autoregressive  model  which  is  then  used  to 
generate  textures. 

5.3  Autoregressive  and  Moving  Average  Models 

In  this  section,  we  will  introduce  the  general  linear 
model  as  defined  by  Box  and  Jenkins  [24J  and  Grabill  [39] 
and  discuss  the  relationships  that  exist  between  it  and 
the  autoregressive  and  moving  average  models.  When  the 
autoregressive  model  is  extended  without  bound  it  is 
essentially  equivalent  to  the  general  linear  model  and  the 
moving  average  model  is  a  subset  of  general  linear  models. 
Allowing  for  a  large  model  reduces  the  complexity  of 
parameter  estimation  and  allows  easy  selection  of  a 
generating  kernel  and  model. 

Many  of  the  equations  in  this  section  also  assume 
that  the  texture  is  one-dimensional  and  that  the 
generating  kernel  is  contiguous.  That  is,  VN+1  follows 


V  .  This  is  consistent  with  notation  presented  earlier  if 
images  are  expressed  in  lexicographic  [40,41]  notation  as 
one-dimensional  vectors.  The  non-cont ig uo us  kernel  may  be 
expressed  as  contiguous  if  the  B^  are  zero  in  the  model. 

The  output  of  a  linear  filter  whose  input  is  white 
noise  may  be  described  using  the  general  linear  time 
series  model 


% 


VN+1  =  eN+l  +  Vn+  *1£N-1+  Vn-2+"* 


=  eN+l+ 


X^jeN-j 


j=l 


(5.5) 


where  {y^+^=VN+^-  y  and  p  is  the  mean  of  the  process, 

assumed  to  be  stationary  here.  Thus  VN+^  is  a  weighted 

sum  of  present  and  past  values  of  the  white  noise  process 

2 

ek  *  ek  usually  has  zero  mean  and  constant  variance  o£. 
The  auto-covariance  is  also  defined  as 


Yt  *  EUkW 


(5.6) 


Notice  that  Eq.  (5.5)  may  be  written  in  a  different 

form 
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(5.7) 


^N+l  *0V  1T1^N-1+  TT2^N-2+'"+  CN+1 


oo 

=  • 
'  1  N-] 


+  e 


N+l 


j=0 


Equation  (5.7)  may  be  derived  from  Eq.  (5.5)  as 


'N+l 


^N+l  ZVN-j 


j=0 


N 


E^j£N-j 


j  =  l 


Therefore 


Vn  =  £n+i+i^o 


V  X^jeN-j 

j  =  l 


+  2-J  jeN-j 
j=l 


eN+l+ti;0VN+  ^  j  eN- j 

j  =  l 


Similarly 


'N-l  ^N-l  S^j£N-j 
j  =  2 


=c>  ^  =  e 

N+l  £N+1 


+*oV(1^0)<'lrN-l+(1'*0)E*jEN-j 

j  =  2 


=  eN+l+*oV,liO)*l 


^N-l 

j  =  2 


+u-*0>£* 

j=2 


(5.8) 


(5.9) 


(5.10) 


(5.11) 


(5.12) 

jEN-j 


106 


oo 


j  =  2 

And  so,  by  continuing  this  process  Eq.  (5.7)  is  found. 

It  may  also  be  shown  that  Eq .  (5.5)  can  be  rewritten 

as 


VN+1  ^(B)eN+l 

where  B  is  the  backward  shift  operator 


(5.13) 


Be.  =  e.  . 
k  k-1 


(5.14) 


B  ek  =  ek-j 


(5.15) 


and 


(MB)  *  (5.16) 

j  =  0 

where  =  1.  ^  (B )  is  often  called  the  transfer  function 

of  the  linear  filter. 

As  certain  constants  or  parameters  must  be  estimated 
from  the  sample  data  available  it  is  sometimes  important 
to  minimize  the  number  of  parameters  required  to 
accurately  represent  a  process.  This  simplifies  analysis 
of  a  model  and  reduces  the  number  of  required 
computations.  The  general  linear  model  containing  an 


107 


injiini  number  of  terms  is  of  little  practical  value. 

Therefore,  if  is  often  expedient  to  allow  the  general 

linear  time  series  process  to  be  reduced  to  a  model  in 

which  the  current  state  of  the  process  may  be  expressed  as 

a  finite  aggregate  of  previous  values  of  the  process  and  a 

driving  error  value  e  .  This  autoregressive  (AR)  model 

N+l 

may  be  written  as 


a. 


K+l  ~  *oV  ^A-l+  ^2^N-2+‘  '  *+  VN-pVl 


(5.17) 


We  may  define  the  auto reg ress ive  operator  <j)(B)  as 


♦  (B)  =  (1-<|)1B-  <P2BZ- 


BP) 


(5.18) 


and  thus  Eg.  (5.17)  may  be  rewritten  as 


VN+1  ^N+l  ’ 


(5.19] 


The  moving  average  (MA)  model  may  be  written  as  a 
special  case  of  Eq.  (5.5)  where  only  the  first  g+1  of  the 
'P  weights  are  non-zero.  This  process  is 

VN+1  =  eN+l  ”  90cN_  0leN-l"  ~  8peN-p 

g  (5.20) 

eN+l“  Z^9jGN-j 
3=0 

As  in  the  case  of  the  autoregressive  model,  we  may  write 
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the  moving  average  operator  as 


0  (B)  = 
and  thus  Eq. 


( 1-0 rtB-0 . B2-  ...  -  0  nBq-0  Bq+1) 
0  1  q-1  q 

(5.20)  may  be  rewritten  as 


8(B)e 


N+l 


(5.21) 


(5.22) 


A  model  incorporating  the  finite-term  concept  of  both 
the  autoregressive  and  moving  average  model  may  be  written 
as 


,u1  =  4>,,v  +<}). .  +  ...+  <J>  +eM  -i 

N+l  ON  1  N-l  p  N-p  N-l 


-  0.£, 0.E.,  i -  •••  ~  0  G.. 

ON  1  N-l  q  N-q 


(5.23) 


or 

^(B)VN+1  =  0 (B) eN+1  .  (5*24) 

This  mixed  autoregressive-moving  average  (ARMA)  processing 
can  be  thought  of  as  the  output  from  a  linear  filter 

whose  transfer  function  is  the  ratio  of  two  polynomials 
0(B)  and  <p  (b  )  when  the  input  is  white  noise  e  . 

An  Example 

* 

To  show  the  relationship  between  an  autoregressive 
process  and  a  moving  process  we  will  consider  the  simple 
example 
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%  _ 

VN+1  ~  eN+l“60eN 


which  is  a  simple  moving  average  process, 
equation  indicates  that 


The 


£N  e0EN-l 


This  implies  that 


CN  =  VeOeN-l 


Substituting  Eq .  (5.27)  into  Eq .  (5.25)  yields 


^N+l  0O^N~eO£N-l+£N+l 


and  by  a  similar  substitution  of 


£  N- 1  ^N-l+0Or'N-2 


into  Eq .  (5.28)  we  see  that 


Vl  =  -‘,O^'0OVN-l-eOrN-2+EN+l 


Continuing  this  process  yields 


N+l 


-E 

j  =  l 


0ov  (N  +  l)  -  j  +  f'N+l 


Thus,  a  finite  moving  average  process  may  be  written 
infinite  autoregressive  process. 


(5.25) 

above 

(5.26) 

(5.27) 

(5.28) 

(5.29) 

(5.30) 

(5.31) 

as  an 
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5.4  Parsimony  Between  Models 


The  above  equations  show  how  the  general  linear  time 
series  model  is  related  to  the  autoregressive,  moving 
average  and  autoregressive-moving  average  models.  The 
concept  of  parsimony  suggests  that  is  is  usually  desirable 
to  express  a  model  with  as  few  terms  as  possible. 
However,  it  has  been  shown  that  i f  we  are  willing  to 
sacrifice  the  concept  of  parsimony  and  deal  with  infinite 
(or  large)  general  linear  models  we  can  use  the  infinite 
autoregressive  model  (Eq.  (5.7)),  also  known  as  the 
general  linear  model,  to  describe  any  general  linear  time 
series  process. 

In  theory,  it  is  possible  to  approximate  as  closely 

•» 

as  desired  any  general  linear  time  series  process  with  a 
finite  a uto reg ress ive  process  of  order  p  by  allowing  p  to 
increase  until  the  desired  closeness  is  obtained.  In 
application,  however,  the  estimation  of  model  parameters 
tk  of  Eq.  (5.7)  may  be  less  accurate  than  desired  due  to 
the  noise  in  the  system.  That  is,  the  noise  of  a 
system  will  cause  error  in  the  estimation  of  the  to  the 
extent  that  some  tt^'s  are  believed  to  be  zero  or  are 
estimated  so  poorly  that  an  inaccurate  model  is  developed. 
It  is  not  clear  to  what  extent  the  accuracy  of  a  model  is 
improved  in  such  a  case  by  using  an  autoregressive  moving 


average  model.  It  is  clear,  however,  that  going  to  such  a 
model  causes  increased  complexity  in  the  parameter 
estimation  process  and  induces  difficulty  into  the 
hypothesis  testing  process  which  is  simple  in  an 
autoregressive  model  case.  Parameter  estimation  for  the 
ARM A  model  requires  a  multiple-pass,  extensive  and 
computational ly-complex  iterative  process  which  is 
complicated  by  the  required  extension  of  our  model  in  a 
two-dimensional  image  to  many  points  and  very  large  sample 
size.  For  this  reason,  the  ARMA  model  was  not  used  for 
texture  simulation  in  this  thesis. 

5.5  Results 

The  linear  (autoregressive)  model  of  Eq.  (5.4)  was 

used  to  simulate  3  variety  of  natural  textures. 

Stationary,  independent  Gaussian  noise  was  used  to  drive 

the  synthesis  process.  The  variance  of  the  noise  was 

estimated  by  applying  the  model  to  the  sample  data  and 

observing  the  prediction  errors.  Images  resulting  from 

the  fitting  of  estimated  models  to  sample  data  are  shown 

later  in  this  chapter.  These  errors,  which  are  often 

called  residuals,  are  pixels  formed  by  tne  difference 

V  -V  ,  where  V  ,  is  the  actual  observed  pixel  value  of 
N+l  N+l  N+l 

tne  sample  parent  image  and  is  the  corresponding 

fitted  value  obtained  by  use  of  the  linear  model.  The 
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standard  deviation  of  these  errors  can  be  measured  and 


1 


used  as  the  standard  deviation  of  pseudo-random 
normally-distributed  noise  in  the  generation  process. 
Actually,  this  information  can  also  be  obtained  during  the 
decomposition  of  the  covariance  matrix. 

The  number  of  pixels  in  each  generation  kernel,  N, 
varied  from  30  to  60.  The  generation  kernels  used  for 
each  simulation  are  shown  in  Table  5.1.  The  simulation 
results  are  shown  in  Fig.  5.1(b)  through  Fig.  5.11(b). 

These  simulations  indicate  that  the  linear  model 
using  stationary  gaussian  noise  produces  acceptable 
simulations  of  a  variety  of  textures  including  grass, 
wool,  leather,  sand  and  water.  As  with  the  binary  model, 
the  simulation  of  bark  (Fig.  5.2)  shows  absence  of  macro 
texture.  The  non-homogeneities  present  in  the  straw 
texture  (Fig.  5.3)  cause  an  "average"  straw  to  be 
generated  which  is  very  similar  to  the  binary  texture 
simulation.  The  cloth  texture  (Fig.  5.4)  is  composed  of 
two  subtextures  and  therefore  simulations  made  using 
statistics  measured  over  the  whole  image  will  be  a  mixture 
of  the  two  subtextures.  The  simulation  of  raffia 
(Fig.  5.11)  has  good  structure  as  the  linear  model  kernel 
is  large  but  sharp  edges  present  in  the  original  texture 
are  absent  in  the  simulation.  The  same  is  true  of  the 
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Table  5.1  TWO-DIMENSIONAL  LINEAR  MODEL  TEXTURE 
GENERATION  KERNELS 
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Figure  5.1  Grass 


115 


(c)  Second-order  Linear 
Model 


(d)  Second-order  Linear 
Model  with  Non¬ 
stationary  Noise 


Figure  5.2  Bark 
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(c)  Second-order  Linear 
Model 


(d)  Second-order  Linear 
Model  with  Non¬ 
stationary  Noise 


Figure  5.5  Wool 


(c)  Second-order  Linear  (d)  Second-order  Linear 

Model  Model  with  Non¬ 

stationary  Noise 


Figure  5.6  Leather 


' 
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(a)  Original  Texture  (b)  First-order  Linear 

Model 


(d)  Second-order  Linear 
Model  with  Non¬ 
stationary  Noise 


(c)  Second-order  Linear 
Model 


Figure  5.7  Sand 


(c)  Second-order  Linear 
Model 


(d)  Second-order  Linear 
Model  with  Non¬ 
stationary  Noise 


(c)  Second-order  Linear 
Model 


(d)  Second-order  Linear 
Model  with  Non¬ 
stationary  Noise 


Figure  5.11  Bubbles 
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simulation  of  sand  when  the  texture  is  examined  in  detail 


on  a  h ig h- reso lut ion  display  device, 


5.6  Rotation  and  Magnification 


The  effects  of  image  texture  rotation  and  scale 
changes  on  the  covariance  function  of  a  texture  can  be 
determined  by  expressing  each  as  a  linear  transformation 
of  the  linear  model.  Using  the  notation  of  the  linear 
model  as  expressed  by  Eq .  (5.1)  we  may  define 


M 


k=l 


(5.32) 


Using  this  notation,  the  maximum  likelihood  estimators  for 
the  mean  and  covariance  matrix  of  N-variate  normal 
distribution,  N(^:u,  T) 


N(*:y,D  =  (2^)'N/2|rr1/2exp[-l/2(^-p)Tr  1(X-u)]  (5.33) 


are 


(5.341 


and 


M 


(5.35) 


k=l 


Relating  this  to  the  facts  discussed  in  sections  3.4,  3.5 


126 


and  ^  ,  we  see  that  the  linear  model  parameter  estimates 
are  obtained  from  the  maximum  likelihood  estimators  for 
the  mean  and  covariance  matrix. 


If  we  then  define  a  linear  transformation  H  such  that 


U) 


k 


(5.36) 


Then 


M 


Ui 


=  -  y 


HX,  =  Hy 


(5.37) 


k=l 


and 


M 


r  = 


-  k  L 


M  ^(HXk  '  HXk>  (HXk  -  HX> 
k=l 

-  s « IX  -  v  <^T  -  *ThT> 


k=l 

M 


(5.38) 


=  'sE'v  V  <xk  -  x,ThT 


k=l 


A  fp 

=  ht  h 
-► 

X 

Thus  if  our  model  is  transformed  linearly,  the  estimates 
required  for  linear  model  parameter  estimation  can  be 
obtained  from  original  model  estimates  and  the 
transformat  ion  itself. 


Rotation  may  be  thought  of  as  a  linear  transformation 
of  coordinate  systems.  Where  our  discrete  image  is 
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considered  to  be  rotated  about  an  axis  by  angle  0  and  the 
standard  row,  column  notation  is  used  (I(n^,n2)  represents 
the  pixel  value  for  image  I  at  row  n^  and  column  ^ ) 


Irotated(nl,n2) 

I  .  .  , (n_sin6+n, cos0 , n, sin6+n„cos0 )  . 

original  2  1  1  2 


(5.39) 


Usually  the  row,  column  addresses  in  the  original  image 
are  fractional.  Therefore,  these  pixel  values  must  be 
specially  defined.  A  widely  accepted  practice  is  to 
estimate  each  value  as  a  function  of  pixels  surrounding 
it.  The  most  likely  candidates  are  nearest-neighbor  and 
bilinear  interpolation [42] .  It  will  be  shown  that  in 
either  case,  the  new  value  may  be  expressed  as  a  linear 
combination  of  values  in  the  original  image. 


A  very  similar  result  may  be  derived  in  scale  changes 
of  a  texture.  Here 


I  i  .  (n. 
scaled  1 


n2)  = 


I  .  .  .  (n.  ■ 

origin? l  1 


a,n_*a) 


(5.40) 


The  origin  of  the  coordinate  system  defines  the  center  of 
the  image  magnification  or  reduction.  In  rotation,  the 
origin  of  the  coordinate  system  defines  the  axis  of 
rotation.  In  both  image  magnification  and  rotation,  the 
row  and  column  addresses  of  the  pixel  in  the  rotated  image 
may  be  fractional.  Again,  the  value  may  be  expressed  as  a 
linear  combination  of  values  in  the  original  image. 
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Consider  the  example  shown  in  Fig.  5.12.  With  angle 
of  rotation  6,  origin  and  magnification  of  4/3,  using 
bilinear  interpolation  (between  rows  the  columns)  we  have 


r 


W1  =  [0.65  V1+0.35  V4] *0.375+[0.65  V2+0.35  Vg]. 0.625 
W_  =  [0.375  V_+0.625  Vc ] • 0 . 35+ [ 0 . 37 5  V-+0.625  V, 1-0.65 

z  zo  ,3  b 

W  =  V 
3  5 

(5.41) 

W4  =  [0.625  V4+0.375  Vy ] • 0 . 65+ [ 0 . 625  V5+0.375  Vg]*0.35 
W5  =  [0.35  V5  +  0.65  Vg]  -0.625+ [0.35  Vg  +  0.65  Vj-0.375  . 


Th  i  s  impl  i  es  that 


(5.42) 
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0.4062b 

0.2167b 

0.0 

0. 2437b 

0. 1312b 

0.0 

0.0 

0.0 

0.0 

0.0 
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0.1312b 

0.0 

0.4062b 

0.2437b 

The  accuracy  of  the  estimate  depends  on  the  ability  oi  the 
interpolating  function  to  accurately  estimate  the  value  of 
off-grid  samples  in  a  discrete  image.  The  covariance 
function  itself  is  being  interpolated  in  this  method. 
Caution  should  be  exercised  when  using  a  nearest-neighbor 
approach  as  the  transfo r.nat  ion  of  a  non-singular 
covariance  matrix  can  be  singular  if  the  rows  of  H  are  not 
independent  vectors. 


By  rotating  covariance  matrices,  we  are  able  lo 
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produce  models  which  can  generate  textures  at  any  angle 
and  magnification  from  one  matrix.  This  could  also  be 
useful  when  trying  to  identify  textures  of  different 
orientation  and  scale  based  on  covariance  statistics. 

5.7  Second-Order  Linear  Model 


When  we 

say  that  a  model  is  a 

linear  or  nonlinear. 

we 

are  referring  to  linearity  or 

nonl inearity 

i  n 

the 

pa  rameters . 

The  value  of  the 

highest  power 

of 

an 

i ndependent 

variable  in  the  model 

is  called  the 

order 

o  f 

the  model. 

For  example. 

Y  =  6lV1+  SnB8+  B0+ 

e 

(5. 

43) 

is  a  second-order  linear  model.  A  general  second-order 
linear  model  with  two  independent  variables  may  be  written 
as 


Y  =e1v1+e2v2+6llv2+622v2+6l2v1v2+60  4-  e  . 


(5.44) 


A  full  second-order  model  with  N  independent  variables 
will  employ  (N*+3N)/2  terms  in  addition  to  the  3 

0 

(constant)  and  e  (error)  terms.  This  general  second-order 
linear  model  may  be  written  as 


Vl  "  6iV62V2  +  '--+8nV60+B11V1  +  612ViV-"*6NNVN+  E 


N  N  N 

=  y^ 6.V.+  y'  y e.  .v.v.+  en+ 

/  j  l  l  L-i  13  l  D  0 

i=l  j=i 


(5.45) 


i=l 
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Second-order  models  have  been  particularly  useful  in 
studies  where  surfaces  must  be  approximated  by  polynomials 
of  low  order.  In  all  cases,  a  second-order  model  will 
"fit"  given  data  as  well  as  or  better  than  a  first-order 
model  that  is  a  subset  of  second-order  models.  This  does 
not  imply  that  the  second-order  model  will  be  more  correct 
however,  as  the  process  which  we  are  attempting  to  model 
may  be  in  fact  a  linear  first-order  process  or  some  other 
type. 


The  use  of  a  second-order  model  to  approximate  the 
surface  of  the  general  stochastic  model  could  have  many 
advantages  over  a  first-order  model.  An  example  of 
fitting  such  a  model  in  one  dimension  to  a  given  set  of 
data  is  shown  in  Fig.  5.13. 

Still  the  linear  first-order  model  may  provide  a  good 


fit  to  the 

data 

and 

the 

mag  ni tud  e 

of  the  unexplained 

variance  in 

the 

data 

may 

be  large 

enough  that  the 

improvement  due  to  the  addition  of  second-order  terms  to 
the  model  may  be  barely  noticeable.  In  two  dimensions, 
the  fitting  problem  is  one  utilizing  a  quadric  surface 
such  as  a  elliptic  paraboloid  or  hyperbolic  paraboloid 
versus  a  plane  to  fit  a  given  set  of  data.  Again,  the  fit 
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may  or  may  not  be  markedly  better.  Adding  second-order 
terms  to  a  model  will  always  produce  a  fit  as  good  as  or 
better  better  than  a  first-order  model  but  the  number  of 
computations  required  to  compute  the  coefficients  and  fit 
the  model  are  much  greater. 

It  is  also  important  to  note  that  the  covariances  of 

the  are  required  in  order  to  obtain  least-square 

estimates  of  the  parameters  3^  in  the  first-order  model 

[20].  Covariance  is  essentially  a  second-order  statistic. 

Therefore,  estimating  the  parameters  of  a  second-order 

model  will  require  the  use  of  fourth-order  statistics. 

Specifically  the  correlation  of  terms  V.  V.  and  V.  is 

needed.  This  may  cause  serious  problems  as  many  cases  the 

variables  in  a  second-order  model  will  be  highly 

2 

intercorrelated .  For  example,  the  terms  V  ,  and 
(if  V  ^  is  highly  related  to  V^)  may  be  strongly 
correlated.  This  situation,  often  referred  to  as 

mul t i-col 1  inea r  i  ty ,  may  cause  problems  during  the 

inversion  or  decomposition  of  the  estimated  correlation 
matrix,  a  necessary  step  in  model  parameter  estimation. 
For  this  reason,  care  should  be  exercised  during  the 
analysis  of  second-order  models. 

Inside  a  circular  radius  of  14  pixels  from  VN+-^  there 
are  307  pixels.  To  search  all  possible  cross  products  in 
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this  region  to  find  the  most  significant  would  require 
over  47,000  cross  products  to  be  examined.  Computation  of 
a  covariance  matrix  containing  all  of  these  terms  is 
impossible  (in  practice).  In  our  study  we  were  limited  to 
investigate  only  820  possible  cross  products  for  entry 
into  the  generation  model.  As  most  of  the  variance  was 
explained  by  the  linear  terms  of  the  model,  most  of  the 
cross  products  were  insignificant  from  a  statistical  point 
of  view.  This  selection  procedure  is  detailed  in  [20]  and 
in  Chapter  3.  Those  that  were  significant  were  entered 
into  the  model  and  a  new  texture  was  generated  using 
Eq.  (5.45)  with  stationary  Gaussian  noise  and  having  zero 
mean  and  fixed  variance  [56]. 

The  results  of  texture  simulations  using  the 
second-order  linear  model  are  shown  in  Figs.  5.1(c)  to 
5.11(c).  On  some  of  these  textures  only  a  slight 
improvement  from  the  addition  of  second-order  terms  may  be 
seen.  In  most  cases,  no  change  can  be  observed  even  when 
the  results  are  displayed  on  a  high-resolution  display 
device.  The  lack  of  improvement  could  be  due  to  the  small 
number  of  cross-terms  examined;  however  we  feel  that  this 
number  is  sufficiently  large  to  show  any  considerable 
improvement  due  to  the  addition  of  second-order  terms  to 
the  linear  model. 
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5.8  Textures  with  Non-Stationary  Noise 


Applying  a  texture  generation  model  to  the  original 
parent  texture  image  data  used  to  estimate  its  parameters 
gives  a  residual  error  image.  When  applying  the  linear 
model  to  a  two-dimensional  texture,  a  two-dimensional 
image  containing  the  pixel  differences  or  residuals 
V„.  -V„.,  is  found.  Here  V,,,.  is  the  prediction  of  the 
next  pixel  in  the  sequence  as  a  linear  function  of  the 
pixel  around  it  according  to  the  model  without  any  noise 
added.  Naturally,  we  would  expect  these  errors  to  be 
small  as  merely  subtracting  one  pixel  from  its 
nearest-neighbor  would  yield  a  small  value  in  most 
natural,  low-noise  images.  Such  an  image  of  residuals  was 
generated  for  the  sand  and  linearly  rescaled  to  show  the 
detail  present  in  the  image  (Fig.  5.14).  Definite 
patterns  are  seen  to  exist  in  this  image  and  thus  a 
violation  of  the  independent  assumption  is  indicated. 
Ideally,  this  residual  image  would  be  uncorrelated  noise. 


A  histogram  showing  the  number  of  VN+^  occurring  at 
each  pixel  value  is  shown  in  Fig.  5.15.  A  plot  indicating 

shown 

in  Fij.  5.16.  As  would  be  expected,  residuals  where  the 


the  mean  of  the  residuals  V„ ,  .  -V,, , ,  versus  V.,,.  i 

N+l  N+l  N+l 


VN  +  1  * s  le;;s  than  0  will  have  a  mean  less  than  zero  and 


those  residuals  where  the  V. 


N+l 


is  greater  than  255  will  be 
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Residual  Sand 


(VN+l"VN+l 


) 


Figure  5.14 


likewise  positive.  Figure  5.17  shows  a  similar  plot  of 

A 

the  standard  deviation  of  the  residuals  versus  V 

N+l 

These  three  figures  seem  to  indicate  that  the  distribution 

/s 

of  the  error  in  the  model  is  related  to  the  value  VN+^ . 
Therefore  the  assumption  of  constant  error  variance  is 


questionable , 


It  may  be  reasonable  to  drive  the 


generation  process  with  noise  which  does  not  have 
stationary  mean  or  variance.  The  effect  of  such  a  change 

in  the  generation  process  was  studied.  Figure  5.18  shows 

/\ 

the  distribution  of  error  VN+1-VN+1  histogram  of  the 
residual  image)  which  appears  to  be  aproximately  normal. 


The  distribution  of  this  error  and  the  relationships 
between  the  predicted  and  actual  pixel  values  was  utilized 
to  generate  textures  using  non-stationary  noise.  The 
procedure  begins  by  generating  a  pixel  VN+1  according  to 
Eq.  (5.45)  excluding  the  error  term.  With  this  predicted 
value  a  random  error  value  e  is  chosen  to  be  added  to 

A 

VN+1*  error  value  e  is  chosen  from  the  distribution 

/v 

of  error  as  a  function  of  V  and  can  have  any  arbitrary 
distribution.  The  next  pixel  will  than  be  computed  in  a 
similar  manner.  Results  of  texture  synthesis  formed  using 
this  model  are  shown  in  Fig.  5.1(d)  through  Fig.  5.11(d). 


The  arbitrary  distribution  of  error  as  a  function  of 
VM4.7  is  calculated  by  applying  the  calculated  linear  model 


Figure  5.17  Standard 

deviation  of 

'»fVi  vs- 


N+l 


Figure  5.18  Histogram  of 

residual  image. 
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to  the  original  parent  texture  and  computing  a  histogram 

/\ 

of  errors  as  a  function  of  ,  . 

N+l 

In  most  cases,  considerable  improvement  is  seen  when 

these  .simulations  are  critically  observed  on  a 

high- resolution  display  device  and  compared  witli  the 

stationary  model  results.  Of  course,  the  information 

required  to  generate  them  is  considerably  greater  also. 

The  distribution  of  errors  as  a  function  of  V  must  be 

N+l 

condensed  and  coded  to  some  degree  to  minimize  storage 

A 

requirements.  For  a  256-grey-level  image  VN+^  usually 
ranges  from  -50  to  305  and  the  errors,  ^n+i-vn+1'  from 
-255  to  +255.  These  ranges  were  determined 

experimentally.  This  would  yield  quite  a  large  amount  of 
data  if  fully  stored.  By  storing  a  small  number  (under 
100),  typical  errors  for  each  range  (and  not  each  single 

A 

value)  of  V  the  number  of  data  values  we  are  required 
N+l 

to  store  can  possibly  be  reduced  to  under  1000.  Therefore 
it  is  believed  that  this  approach  of  using  non-stationary , 
non-Gaussian  noise  to  generate  textures  may  be  quite 
acceptable  even  with  severe  storage  limitations. 

5.9  Conclusion 

The  results  in  this  chapter  indicate  that  many 
natural  textures  are  well  simulated  using  a  large 


autoregressive  model.  Adding  second-order  terms  to  the 


model  improves  the  results  slightly  but  the  resulting 
increase  in  computational  complexity  makes  this 
second-order  simulation  method  difficult  to  implement. 
Using  non-stationary  noise  in  the  generation  process 
improves  the  simulations  considerably  when  the  textures 
are  viewed  critically.  The  subsequent  increase  in  storage 
and  computation  required  by  this  addition  is  small  and 
therefore  this  model  modification  should  be  considered  in 
most  applications. 
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CHAPTER  6 


MULTIPLE  MODEL  TEXTURE  GENERATION 


G.l  Introduction 


In  this 

chapter , 

we  will  present  three 

methods 

o  f 

generating 

textures 

using  multiple  texture 

models . 

The 

first  method 

introduces  a  set  of  generation 

kernels 

that 

is  used  to  synthesize  texture  pixels  in  a  multiple-pass 
manner.  Associated  with  each  of  these  kernels  is  a  unique 
model.  This  method  could  be  useful  in  generating  textures 
which  have  very  coarse  structure.  The  second  method  uses 
a  piecewise-1 inear  method  of  fitting  the  model  to  parent 
texture  data.  The  model  chosen  during  the  generation 
process  is  allowed  to  depend  on  the  pixel  values  in  the 
kernel.  Although  the  fit  of  the  model  is  better,  the 
synthesis  results  show  little  improvement  over  the 
single-model  approach  with  a  linear  model.  The  third 
method  of  texture  generation  presented  in  this  chapter 
uses  a  additional  image  to  determine  the  model  number  to 
be  used  during  the  generation  process.  This  composite 
generation  method  could  be  useful  when  a  texture  is 
actually  composed  of  a  set  of  subtextures  as  it  allows  a 


unique  model  for  each  of  these  subtextures  to  be  used  in 
the  generation  process. 


6.2  Sk ip-Generate  Method 

Simulating  textures  which  have  a  fine  structure  is 
usually  a  much  easier  process  than  simulating  textures 
with  coarse  structure.  This  occurs  because  the  linear 
model  contains  fewer  terms  if  the  texture  pixels  become 
uncorrelated  over  a  small  distance.  For  the  same  texture 
at  a  greater  magnification,  the  pixels  become  highly 
correlated  and  the  linear  model  will  be  forced  to  contain 
more  terms.  As  the  texture  becomes  more  coarse,  more 
time-consuming  statistical  measurements  must  be  taken  on 
the  parent  texture  over  larger  windows.  Motivated  by 
these  problems,  the  texture  generation  algorithms  in  this 
section  have  been  developed. 

In  the  texture  work  so  far,  pixel  V  was  generated 
based  on  pixels  above  or  to  the  left  of  it  (see 
Fig.  3.1(b)).  As  discussed  in  Chapter  3,  the  kernel  does 
not  have  to  be  contiguous.  This  kernel  shape  is  chosen  to 
insure  that  the  image  space  of  our  synthesized  texture  was 
filled  during  the  generation  process.  However,  generating 
pixels  along  a  row,  row  by  row  is  not  the  only  way  of 
filling  an  image  space. 
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Consider  the  non-contiguous  kernel  mask  in  Fig.  6.1. 
If  the  spacing  between  the  pixels  in  this  mask  is  8,  using 
the  linear  model  in  Eq .  (5.4)  to  generate  the  right-most 
pixel  in  the  bottom  row,  we  can  generate  every  8th  pixel 
along  every  8th  row.  At  each  step  the  next  pixel  is 
generated  based  on  the  previously-generated  pixels  around 
it  (ignoring  boundary  conditions).  After  generating  an 
image  with  this  type  of  spacing,  the  pixels  midway  between 
the  previously-generated  pixels  on  each  row  may  be 
generated  using  the  mask  in  Fig.  6.2.  In  this  mask,  the 
pixel  with  the  "x"  in  it  denotes  the  next  pixel,  V  , ,  to 
be  generated  according  to  Eg.  (5.4).  Naturally,  the 
linear  model  used  in  this  step  will  have  different 
coefficients  than  the  previous  one.  It  is  also 
interesting  to  note  that  new  pixels  depend  not  only  on 
previously  generated  pixels  above  them  (as  with  the  mask 
in  Fig.  3.1(b))  but  depend  also  on  the  pixels  below  them. 
Still,  ignoring  boundary  conditions,  each  pixel  depends 
only  on  previously  generated  pixels.  At  the  next  step  a 
mask  similar  to  that  in  Fig.  6.3  can  be  used  to  fill  in 
the  pixels  midway  between  the  previously-generated  pixels 
in  each  column.  Again  pixels  are  allowed  to  depend  on 
pixels  around  them. 

By  repeatedly  using  the  masks  in  Fig.  6.2  and 
Fig.  6.3  with  successively  closer  and  closer  pixel 


143 


□ 


□ 


□ 


□ 


□ 


□  □  □  □  □ 

□  □  □  □  □ 

□  ns 

Figure  6.1  First-pass  Generation  Kernel 
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Figure  6.2  Second-pass  Generation  Kernel 
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spacing,  the  texture  simulation  image  space  is  filled.  An 
example  showing  the  pixels  generated  at  each  successive 
pass  is  shown  in  Fig.  6.4.  Wore  importantly,  to  determine 
the  linear  model  for  each  mask,  only  one  covariance  matrix 
is  required  and  can  contain  as  many  or  as  few  terms  as 
desired.  The  process  of  collecting  statistics  for  one 
matrix  is  not  beyond  the  complexity  that  we  would  want  to 
undertake  for  the  small  number  of  times  required  by  this 
process.  Naturally,  any  other  stochastic  process  may  be 
substituted  for  the  linear  model.  As  before,  only  the 
measurements  required  to  estimate  the  parameters 
corresponding  to  each  mask  need  to  be  taken.  This  number 
depends  on  the  spacing  of  the  pixels  in  the  first  mask, 
which  should  be  a  power  of  two.  Other  odd-shaped  kernels 
and  kernels  whose  spacing  is  not  a  power  of  two  could  be 
designed  to  form  space- f i 1 1  ing  sets.  Most  would  require 
more  models  to  be  estimated  and  would  provide  little 
additional  information. 

Texture  simulations  using  this  method  are  shown  in 
Figs.  6.5-6.12.  Only  a  slight  improvement  is  seen  in  some 
of  the  texture  simulations  over  the  synthesis  done  by  the 
earlier  single  linear  model.  Most  of  these  textures  are 
apparently  well  simulated  by  a  carefully  chosen  model  and 
the  results  are  not  critically  dependent  on  the  coarseness 
of  the  textures. 
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Figure  6.3  Third-pass  Generation  Kernel 
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Figure  6.4  Filled  Space  of  Skip-generate  Method 
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Figure  6.9  Skip-generate 
Leather 


Figure  6.10  Skip-generate 
Sand 


Figure  6.11  Skip-generate 
Water 


Figure  6.12  Skip-generate 
Wood 
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A  word  of  caution  should  be  added  concerning  the 
computations  involved  in  the  linear  model  coefficient 
calculation  of  this  method.  During  the  later  stages  of 
the  skip-generate  method,  the  pixels  in  the  generation 
kernel  become  highly  correlated  as  the  distance  between 
them  decreases  with  each  pass.  This  may  cause  the 
correlation  or  covariance  matrix  of  the  model  to  be 
ill-conditioned.  To  avoid  numerical  problems,  the  number 
of  variables  entered  into  the  process,  and  therefore  the 
number  of  steps  involved  in  the  matrix  decomposition 
process,  should  be  kept  to  a  minimum  in  some  cases.  The 
use  of  ridge  regression  techniques  [43,44,46]  might  also 
be  considered. 

6.3  Piecewise-Linear  Models 

When  generating  textures  using  the  general  linear 

model  described  by  Eq.  (5.4)  and  the  generating  kernel  in 

Fig.  3.1(b)  the  same  model  is  used  regardless  of  the 

values  of  the  pixels  V^,...,VN»  By  developing  more  than 

one  linear  model  and  allowing  the  choice  of  the  model  at 

each  pixel  generation  step  in  the  synthesis  process  to  be 

dependent  on  some  functional  value  of  V, ,...,V„, 

1  N 

F(V,..., )  a  new  synthesis  model  is  formed. 

To  illustrate  this  concept  consider  the  data  in 
Fig.  6.13(a).  If  we  were  to  fit  one  linear  model  to  the 
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data  in  order  to  predict  from  it  would  look  like  the 
single  line  running  through  the  data  in  Fig.  6.13(a). 
This  linear  model  could  then  be  used  to  predict  V 2  based 
on  the  value  of  V  .  But  if  we  allow  the  choice  of  our 
linear  model  to  be  dependent  on  the  value  of  ,  then  for 
an  incoming  value  of  we  choose  a  model  whose  domain 
includes  V-^  to  predict  V2.  For  5  linear  models,  the 
straight  lines  are  shown  in  Fig.  6.13(b) .  The  fit  to  the 
data  using  multiple  linear  models  will  always  be  as  good 
as  or  better  than  that  of  the  single  linear  model.  That 
is,  the  mean  square  error  will  generally  be  reduced  using 
multiple  models. 

Using  multiple  linear  models  for  texture  synthesis  we 

would  generate  pixels  V  based  on  pixels  V  ,...V  in  the 

N+l  1  N 

following  way.  First,  we  compute  a  function,  F,  of  the 
V^,...,V  pixels  which  allows  us  to  choose  the  proper 
linear  model.  Then  using  this  model  with  the  values 
V^,...,VN  we  predict  VN+^  and  add  noise.  This  process  is 
diagramed  in  Fig.  6.14. 

Ideally,  the  function  F  should  be  chosen  to  minimize 
the  total  mean  square  error  resulting  from  fitting  the 
limited  number  of  models  to  the  sample  data.  This  is  very 
difficult  to  do  in  practice  however  as  for  M  larger  than  3 
we  are  fitting  multiple  hyperplanes  to  data  in  an  N+l 
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dimensional  space. 


One  texture  synthesis  of  sand  was  done  using  the 
multiple  linear  model  (see  Fig.  6.17).  In  this  case  eight 
models  were  used  and  the  model  number  was  chosen  by 
examining  the  pixel  immediately  to  the  left  of  the  pixel 
being  generated.  The  range  of  this  pixel,  0  to  266,  was 
divided  into  8  equal  subranges  and  the  model  was  chosen 
according  to  the  subrange  into  which  the  value  fell.  Only 
a  slight  improvement  over  the  single  linear  model 
synthesis  (see  Chapter  5)  is  seen  even  though  the  same 
kernel  shape  was  used.  No  other  simulations  have  been  run 
using  this  model  as  it  is  felt  that  little  improvement 
will  result.  Also  the  linear  model  containing 
cross-product  terms  in  Chapter  5  probably  provides  a  very 
good  fit  in  most  cases  and  in  more  dimensions.  In  one 
dimension  the  model  of  Chapter  6  would  fit  the  data  in 
Fig.  6.13  with  a  quadratic  curve. 

6.4  Field-Definition  Stochastic  Model 

Another  method  of  using  multiple  stochastic  models  is 
to  generate  an  image  of  fields  defining  the  model  number 
to  be  used  in  a  second  pass.  Such  an  approach  would  be 
useful  in  simulating  textures  which  have  multiple 
sub-textures  within  them.  A  simple  analytical  example  is 
shown  in  Fig.  6.15.  Real  world  examples  might  include 


such  things  as  a  brick  wall  where  the  texture  of  the 

bricks  is  different  than  the  texture  of  the  concrete 
separating  them.  It  was  felt  that  this  type  of  approach 
might  be  useful  in  the  simulation  of  bark  (see 
Fig.  3.6(a))  which  has  a  strong  macro  structure.  A  method 
to  separate  this  texture  into  two  fields,  which  would 
later  define  the  model  to  be  used,  was  designed.  This 

result  (Fig.  6.18)  was  obtained  by  successively  passing 
smart  median  filters  of  varying  sizes  over  a  binary  image 
(which  was  obtained  by  thresholding  an  original  continuous 
grey-level  image)  (see  Fig.  5.2(a)). 

The  smart  median  filters  replace  the  center  pixel  of 

a  window  with  the  median  only  if  certain  conditions  are 

net.  The  window  is  passed  over  the  entire  image  pixel  by 

pixel  along  a  row,  row  by  row  in  a  two-dimensional 

convolution  manner  as  in  Fig.  6.15.  Let  I  (n  ,n  )  be  the 

112 

input  image  and  I  (n  ,n  )  be  the  output.  Let  the  pixel 

0 

values  be  0  (black)  and  1  (white) .  Let  N  denote  the 

B 

number  of  black  pixels  in  the  window  being  processed  with 
center  anc^  let  be  the  number  of  white  pixels 

in  the  window.  For  the  binary  case, “the  smart  median 
filter  is  defined  as 
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Figure  6.15  Analytic 

Sub-textures 


Figure  6.16  Two-dimensional 
Convolution 


r 


0, 

if 

nb 

VNw 

N 

> 

Thresh 

1, 

if 

W 

VNw 

> 

Thresh 

(n^jn^)/  otherwise 

(6.1) 


Equation  (6.1)  indicates  that  the  center  pixel  of  the 
window  is  replaced  by  the  median  if  the  percentage  of 
black  or  white  pixels  is  above  a  specified  threshold.  The 
term  "replaced"  is  used  loosely  in  this  context  to 
indicate  the  visual  appearance  of  replacement  when  the 
input  and  output  images  are  superimposed.  In  the  binary 
case,  the  median  of  the  pixels  in  a  window  is  equal  to  the 
value  of  the  most  frequently  occurring  pixel  in  the 
window.  The  threshold  and  window  size  varied  in  the 
successive  passes  over  the  image.  The  window  sizes  and 
thresholds  for  each  of  the  passes  used  to  obtain  Fig.  6.16 
were  1-0.50,  3-0.50,  5-0.50,  7-0.75,  11-0.78,  7-0.50. 
This  mul t i pi e- pass  procedure  helped  to  retain  detail  and 
eliminate  fields  too  small  for  useful  measurements. 


The  field-definition  stochastic  model  was  also 
applied  to  the  non- sta t iona r y  cloth  texture  (see 
Fig.  5.14(a)).  This  texture  is  clearly  composed  of  two 
alternating  subtextures.  In  this  case,  the  texture  fields 
can  be  extracted  by  hand  (see  Fig.  6.21). 
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Figure  6.19  Field-definition  Figure  6.20  Field-def ini- 
for  Bark  tion  Bark 

Simulation 


Once  the  two-field  images  for  bark  and  cloth  were 
obtained  they  were  processed  again  to  define  a  total  of 
four  fields,  the  two  original  fields  and  two  border  or 
transition  fields.  The  transition  fields  are  required  as 
the  kernel  of  points  over  which  relational  measurements 
(usually  second-order  statistics)  are  taken  may  not  fall 
wholly  within  one  of  the  two  fields.  Such  border 
measurements  will  surely  increase  the  inaccuracy  of  the 
models  for  the  original  fields.  So  these  border  regions 
are  considered  to  be  unique  texture  fields.  They  are 
mathematically  defined  by  passing  the  original  linear 
model  texture  generation  kernel  obtained  in  Chapter  5  over 
the  texture  and  defining  the  relative  importance  of  each 
pixel  to  be  equal  to  the  absolute  value  of  its  associated 
6^  coefficient  of  the  linear  model  (Eq.  (5.4)).  These 
coefficients  are  usually  largest  near  the  eye  of  the 
kernel.  The  "importance"  of  the  eye  itself  was  set  to 
1.5  max(B^).  If  a  large  percentage  (90%  for  bark,  95%  for 
cloth)  of  the  total  "importance"  fell  within  one  of  the 
two  original  fields,  the  point  would  be  considered  as  a 
member  of  that  field.  If  not,  the  point  would  become  a 
member  of  the  transition  field  based  on  its  position  with 
respect  to  the  two  fields,  A  and  B.  If,  at  the  point 
being  processed,  A  is  to  the  left  of  B  then  the  point  is 
in  one  transition  field,  if  B  is  to  the  left  of  A  then  the 
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point  is  in  the  other  transitional  field.  The  right-left 
orientation  is  used  as  the  texture  generation  process  is 
usually  done  in  a  1 ef t- to- r ig ht  manner  along  each  image 
row. 


The  results  from  processing  the  original  field  data 
for  the  textures  bark  and  cloth  are  shown  in  Fig.  6.19  and 
Fig.  6.22.  The  four  fields  are  represented  using  four 
grey-levels.  The  black  regions  on  the  border  may  be 
considered  to  be  undefined  as  the  kernel  points  do  not  lie 
with  the  image  region  in  these  areas. 

Tnese  field  definitions  were  used  both  in  the 
statistics  gathering  process  as  well  as  the  texture 
synthesis  process  using  the  four  calculated  models.  The 
models  were  linear  (see  Eq .  (5.4))  and  the  methods  used  in 
their  estimation  were  discussed  earlier  in  Chapter  5.  The 
synthesis  results  are  shown  in  Fig.  6.20  and  Fig.  6.23. 
Unfortunately,  the  field  structure  is  not  strongly 
apparent.  The  primary  reason  is  that  each  model  requires 
a  number  of  pixel  generations  before  it  reaches  a  steady 
state  and  in  most  cases  this  is  much  greater  than  the  size 
of  most  of  the  fields  in  the  bark  texture.  By  observing 
the  border  regions  of  many  synthesis  runs,  it  seems  that  a 
steady  state  is  reached  after  10-20  pixels  have  been 
previously  generated.  In  the  cloth  texture  synthesis  of 
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Fig.  6.23,  the  randomness  of  the  noise  seems  to  destroy 
the  structure  of  the  texture  quite  frequently.  This  is 
illustrated  in  Fig.  6.24  where  only  one  of  the  four  models 
was  used  to  generate  the  entire  texture.  Each  of  the 
models  contained  the  same  set  of  kernel  points  used  in  the 
original  linear  model  of  Chapter  5.  It  is  possible  that 
improvement  could  be  made  by  choosing  a  different  set  of 
points  for  each  field-model. 

In  an  actual  simulation,  once  a  field  image  is 
obtained  a  method  must  be  developed  to  simulate  this  field 
texture  and  this  field  texture  will  then  in  turn  be  used 
to  choose  the  model  numbers  in  the  generation  of  final 
synthesis.  Generating  textures  with  only  a  few  grey 
levels  can  be  done  using  more  complete  stochastic 
statistics,  perhaps  N-grams,  but  the  large  size  of  the 
fields  may  require  that  a  method  such  as  the  ski p-generate 
method  (discussed  in  section  6.2)  be  used.  For  many 
textures,  such  as  cloth,  the  field  generation  process 
would  be  simple. 

More  work  should  be  done  in  the  multiple  model  area 
as  the  storage  requirements  for  such  models  is  very  small 
but  can  potentially  produce  improved  results.  A  great 
number  of  combinations  and  approaches  are  possible  in  this 
area  and,  in  many  cases,  may  be  chosen  by  the  application 
or  the  particular  textures  being  simulated. 
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CHAPTER  7 


8EST-F IT  TEXTURE  MODEL 


7 . 1  I n t  rod  uc  t i on 

A  method  of  generating  texture  simulations  according 
to  their  Nth  order  densities  was  investigated  for  binary 
textures  in  Chapter  3.  The  simulations  resulting  from 
this  Markov  process  resembled  their  parental  textures 
quite  closely  in  most  cases.  When  applying  a  similar 
concept  to  multi-grey  level  imagery,  the  limits  of 
computer  storage  are  soon  reached.  To  circumvent  this 
constraint,  a  new  method  of  texture  synthesis  was 
developed  and  applied  to  a  number  of  textures.  Simulation 
results  using  this  method  are  given  in  this  chapter. 

7.2  N-grams  in  Continuous  Imagery 

In  binary  texture  generation  based  on  N-grams  a 
single  functional  value  P  (V  , /V,  ,  .  .  .  ,  V„  )  was  stored  for 
each  possible  pattern  (V^  ,  ,  .  .  .  ,  )  where  the  \A  '  s  can  be 

zero  or  one.  This  value,  also  called  a  generation 
parameter,  represented  the  conditional  probability  that 
the  next  pixel,  ,  in  the  generation  process  would  be  a 


zero-valued,  black  pixel.  The  '  s  were  chosen  by  a  best 

linear  model  fit  detailed  in  Chapter  3  and  therefore  the 

kernel  of  previous  pixels  (V^,...,V  )  is  not  necessarily 

contiguous  (see  Figure  7.1).  Details  concerning  the 

estimation  of  P  (Vj^+l /V1  '  *  *  '  '  VN  )  ^rorn  a  parent  texture  are 

given  in  Chapter  3.  For  binary  textures,  this  single 

value  is  sufficient  to  define  the  distribution  of  . 

N+l 

given  V  , ...,V  .  The  number  of  different  functions  which 
1  N 

N 

must  be  stored  is  2  .  In  the  generation  process  each 
pixel  is  generated  based  on  the  values  of  the  pixels 

V  , . . . ,VN  surrounding  it  and  on  a  computer-generated 
uniformly-distributed  random  variable.  The  texture 
simulations  are  generated  pixel  by  pixel  along  a  row  until 
each  row  is  complete.  Pixel  generation  along  the  edges  of 
an  image  can  be  handled  in  a  variety  of  ways  but  in 
Chapter  3  pixels  in  these  border  regions  were  assumed  to 
be  any  random  value,  0  or  1,  if  they  were  outside  the 
image  boundaries. 

A  similar  approach  could  be  used  to  generate 

multi-grey-level  textures.  For  a  texture  containing  g 
N+l 

grey  levels,  g  different  functions,  P (VN+^/V^, . . . , VN) , 

roust  be  stored.  (Actually  only  (g-1)  gN  are  required  as 
g-1 

P  (X/V^,  .  .  .  ,V^) -1  for  all  V^).  Storage  limitations 

are  soon  reached.  Also  estimation  of  P(V„  , /V V  )  is 

N+ 11  N 

difficult  as  multiple  occurrences  of  the  pixel  pattern 
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V  / • • • / V  may  not  exist  in  the  parent  texture.  Therefore 
1  N 

even  without  storage  limitations  tne  problems  of 

estimating  P  (V  /V  , .  ..,V  )  from  a  given  parent  texture, 
N+l  1  N 

which  represents  the  distribution  of  V  given  the  values 

N+l 

of  V  , . . . , V  is  complex. 

1  N 

This  estimation  problem  no  doubt  has  a  number  of  ad 

hoc  solutions.  The  problem  is  basically  that  for  large  N 

and/or  large  g,  there  may  not  be  a  suitable  number  of 

occurrences  of  the  pattern  V^,...fV  to  adequately 

estimate  the  distribution  P  (V  /V  ,...,V  )  given  a  finite 

N+l  1  N 

sample  size.  Even  though  a  certain  pattern  never  occurs 

or  rarely  occurs  in  our  sample  parent  texture  it  is  not 

implied  that  such  a  pattern  is  impossible  and  will  never 

occur  in  our  simulation  synthesis.  We  might  often  find 

numerous  occurrences  of  this  pattern  if  our  sample  size  or 

the  size  of  our  parent  texture  was  increased,  especially 

in  noisy  and  fine-structured  textures.  But  as  this  very 

large  sample  may  not  be  present,  we  must  estimate 

P (V  /V  ,...,V  )  for  all  V  ,...,V  based  on  available 
N+l  IN  IN 

samples. 

One  approach  would  be  to  use  sample  patterns  which 

closely  resemble  but  which  may  not  be  exactly  the  same  as 

each  pattern  (V  ,...,V  ).  That  is  in  a  pictorial  sense, 
1  N 

we  use  patterns  of  (V_,...,V  )  which  look  "close  to"  the 

1  N 


pattern  for  which  we  are  attempting  to  estimate 

P (V  /V  f • • - #V^) .  Therefore  samples  in  our  sample  parent 

texture  may  be  used  to  estimate  numerous  P  (V., , , /V,  . . . .  ,V„) 

N+l  1  N 

and  not  just  those  they  fit  exactly.  The  concept  of  a 

distance  function  must  be  used  to  numerically  define 

"close  to".  Given  two  patterns,  one  from  our  sample 

texture  and  the  other  from  the  conditional  probability  of 

the  kernel  we  are  attempting  to  estimate,  the  distance 

measure  can  be  used  to  determine  the  value  of  that  sample 

in  estimating  P (V  /V  ,...,V  ).  If  the  fit  between  the 

N+l  1  N 

kernel  pattern  and  the  pattern  in  the  sample  texture  is 

good  the  associated  value  of  V^+1  in  the  parent  texture 

will  be  valuable  in  estimating  P (V  /V  ,...,V  ). 

N+l  1  N 

Normally,  when  N  and  g  are  small  or  when  we  have  many 
samples  for  any  given  we  can  use  the  histogram 

of  the  associated  VN+^  to  estimate  P  (VN+^/Vj , . . . , VN) . 
Here  the  relative  number  of  times  a  particular  value  of 
occurs  given  a  pattern  indicates  the  conditional 
probability  we  are  attempting  to  estimate.  This  was 

discussed  in  section  3.2.  tvnere  a  distance  measure  is 
used  instead,  a  good  fit  could  be  considered  to  be 
synonymous  wi th  high  frequency  of  occurrence  of  that 
pattern  and  a  poor  fit  with  low  frequency  of  occurrence. 

If  such  a  method  of  estimating  these  conditional 
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probabilities  is  used  we  are  still  faced  with  a  huge 
storage  problem.  For  this  method  to  be  practical,  the 
storage  requirement  n.ust  be  reduced.  From  an  information 
standpoint,  it  is  interesting  to  note  that  a  method  of 
estimating  N-grams  or  conditional  probabilities 
P (VN+i/Vi , . . . ,VN)  from  a  sample  parent  texture  image 
produces  gN+^  data  values  from  M  pixel  samples  where  M  is 
the  size  of  the  square  parent  texture  image  in  pixels. 
For  large  g  and  N  this  is  a  drastic  increase  in  data.  But 
the  actual  information  content  can  really  never  be  greater 
than  that  content  of  the  sample  parent  texture  image. 
Therefore,  this  M  value  represents  an  upper  bound  on  the 
amount  of  data  we  should  use  to  generate  a  texture 
simulation.  Any  amount  of  data  exceeding  this  will 
contain  redundant  data. 

7.3  The  Synthesis  Method 

Combining  this  concept  of  upper  bound  with  the  idea 
of  forming  a  distance  measure  to  compare  two  texture 
kernel  patterns  leads  to  a  new  texture  synthesis  method. 
In  this  method,  we  generate  the  next  pixel  based  on  the 
pixels  in  the  kernel  surrounding  it  (see  Figure  7.1)  and 
their  comparison  to  similar  kernels  in  the  parent  texture. 
This  comparison  is  made  by  passing  the  kernel  currently 
present  in  the  simulation  process  over  the  parent  texture 
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and  computing  the  distance  function  at  all  possible  points 
(see  Figure  7.2).  Denoting  the  pixels  in  the  parent 
texture  by  X  ,  i  ,  j=0  ,  .  .  .  , /M-l  and  the  pixels  in  the 

i/ j 

kernel  V.,...,V  by  Y.  .  ,  we  can  compute  a  comparison 

-L  N  1  r  J 

image 

c  ,  =  COMPARISON  (X .  ,  .  ,Y.  .)  n  ,, 

a,b  i+a,g+b'  i,]  (7.1) 

for  all  a,b  such  that  the  kernel  is  within  the  boundaries 
of  the  texture. 

One  possible  comparison  function  would  be 

correlation.  Assuming,  without  loss  of  generality,  that 

our  kernel  is  contiguous  as  in  Figure  7.3  and  the  elements 

are  denoted  as  Y.  .,  this  function  would  be  defined  as 
1  >  3 


The  problem  with  this  particular  distance  measure  is  quite 
serious.  Correlation  does  not  take  into  account 
differences  in  over-all  mean.  For  example,  the  kernels  in 
Figure  7.3  are  perfectly  correlated  but  their  means  differ 
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significantly.  Differences  in  first-order  statistics 
between  these  kernel  patterns  will  not  be  detected  by  a 
correlation  measure  and  so  another  comparison  function  to 
supplement  a  correlation  value  would  be  required. 

A  better  comparison  function  would  be  the  mean  square 
difference  (MSD) .  This  is  defined  as 

MSD  ,  =  /  /  (X.^  .  ,,-Y.  ,)2 

a,b  L-ul—i  i+a,g+b  i,j  (7.3) 

i  j 

where  i,j  must  be  within  the  coordinate  range  of  the 
kernel  as  in  Eq .  (7.1).  This  comparison  function  will 
detect  differences  in  first-order  statistics  between  two 
pixel  patterns  (such  as  those  in  Figure  7.3)  as  the  MSD 
function  is  a  sum  of  squares  of  differences.  Whereas  the 
correlation  coefficient  of  Eq .  (7.2)  varies  between  -1.0 

and  1.0  and  is  largest  when  the  fit  is  good,  tne  MSD 
measure  is  small  when  the  fit  is  good  and  it  is  always 
positive. 

The  MSD  function  weights  the  comparison  of  all 
elements  in  a  kernel  equally.  Having  studied  many  texture 
generation  models  we  immediately  recognize  that  this  fit 
is  not  properly  weighted.  The  few  pixels  which  are 
closest  to  Ynext  in  proximity  are  far  more  important  when 
predicting  ynext  t^ian  those  which  are  far  away.  So 


Eq.  (7.2)  must  be  modified  to  give  the  weighted 

mean-square  difference  (WM5D) 

WMSD  ,  =  y"'^T'(X._1_  ...-Y .  .  )2.  W.  .  .  (7.4) 

a,b  i+a,  j+b  i,j 

i  j 

A  possible  choice  for  W  is 


W.  . 
i»D 


(7.5) 


2  2 
^i_iNEXT^  +^“^NEXT^ 


where  R  is  the  euclidean  distance  between  pixel  Y.  .  and 

1/3 

the  kernel  eye  and  the  coordinates  of  the  eye  are 

NEXT  1 


g iven  by  ‘ t NEXT - 3 NEXT> - 


As  the  first  step  in  comparing  a  given  kernel  Yi  ^  to 

all  kernels  in  the  parent  texture,  for  each  point  (a,b)  in 

the  parent  texture,  ignoring  edges,  the  WMSD  is  computed 

resulting  in  an  image  of  WMSD's.  Where  the  fit  between 

the  generated  kernel  Y.  .  and  the  image  X.  .is  good,  we 

1  /  J  1  f  D 

would  expect  WMSD^  ^  to  be  small.  The  smallest  WMSD 
represents  the  "best"  fit  according  to  our  norm.  We  could 
choose  the  ^jjEXT  associated  with  this  best  fit  at  point 
(a,b)  to  be  our  next  pixel  in  the  generation  process, 
however  this  can  cause  problems.  First  of  all,  the 
generation  process  would  "lock  in"  on  the  parent  texture 
and  the  generated  texture  could  very  well  become  just  an 
exact  copy  of  the  input  parent  texture.  Second,  we  know 
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p 


ideally  that  Y  has  a  distribution,  not  just  a  mean. 

NEXT 

In  the  autoregressive  model  of  Chapter  5  we  gave  YNEXT  a 

distribution  by  adding  random  noise  to  it.  Although  this 

could  be  done  here,  such  an  approach  would  fail  to  use 

additional  information  contained  in  the  WMSD  image.  There 

may  be  a  set  of  points  (a,b),  all  exhibiting  a  good  fit  to 

the  kernel  pattern  Y.  ..  In  fact,  the  best  fit  may  have  a 

i/3 

noisy  Y  and  the  other  good  fits  could  provide 

NEXT  1 

information  to  improve  the  prediction  of  the  YNEXT  in  the 
generation  process.  Using  a  set  of  best  fits  is 
equivalent  to  increasing  our  sample  size.  We  look  at  a 
set  of  similar  patterns  to  pick  our  Y 

NEXT 


At 

this  point  there  are 

numerous 

ways 

to  proceed. 

Logically  those  patterns 

with  the 

"  best 

"  fit  should 

prov ide 

better  estimators 

for  yNEXT 

so 

some  kind  of 

weighting  decision  is  needed  to  choose  the  relative 

importance  of  the  WMSD's  found.  If  we  search  through  the 

WMSD  image  and  find  the  minimum  value,  WMSD  .  ,  and  scale 

min 

all  the  WMSD's  by  that  we  form  a  new  image  MAXI 


MAXI 


a,  b 


WMSD  . 
min 


WMSD 


a,b 


(7.6) 


This  image  has  the  value  1.0  at  the  best  fit  point  and 
values  0  <  MAXI  <  1.0  elsewhere. 
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Here  we  can  look  at  the  MAXI (a, b)  image  and  study  its 

range.  I£  0.16  MAXI  <  1.0  it  is  implied  that  the  worst 

fit  yields  a  0.16  value.  Somehow  that  worst  fit  should  be 

translated  to  imply  that  the  probability  of  choosing  the 

YNEXT  associated  with  that  point  (a,b)W0RST  *s  near^y  0.0. 

The  simplest  way  of  doing  that  is  to  take  powers  of  the 

image  MAXI  (a, b).  The  maximum  remains  1.0  while  smaller 

numbers  approach  0.0.  For  example  (1. 0)^=1. 0  but 
10  -8 

(0.16)  =1x10  .  We  do  this  to  obtain  an  jad  hoc  estimate 

of  P^YNEXT//Yi  j)*  After  experimentally  studying  the 
values  of  MAXI  (a, b)  and  its  powers,  the  value  of  16  was 
chosen  and  a  new  image  PDFUNS 

PDFUNS  ,  =  (MAXI  .  )16 

a ,  b  a,b  (7.7) 

was  used  to  estimate  the  probability  density  function 
P(^Ht.vm/Y;  -;)•  The  values  in  the  PDFUNS  image  are 

Nila  I  1 1  J 

generally  very  small  with  less  than  1%  of  the  image  points 
having  value  greater  than  0.1.  As  a  rule  of  thumb,  it  can 
be  argued  that  1%  to  0.05%  of  the  128x128  PDFUNS  values 
should  be  greater  than  0.1.  A  larger  percentage  would 
increase  undesired  randomness  and  noise  in  the  synthesized 
image  and  a  smaller  number  could  cause  "lock  in"  on  the 
parent  texture.  The  value  16  was  also  chosen  for 
convenience  and  computational  efficiency  as  it  can  be 
computed  with  only  4  multiplications  and  minimal  data 
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storage.  More  study  could  be  made  concerning  the  effect 
of  the  value  of  this  variable  on  the  simulation  results 
and  other  approaches  to  creating  a  PDFUNS  image  from  the 
MAXI  values  could  be  tried.  Studying  a  histogram  of  MAXI 
values  might  also  be  very  informative. 


PDFUNS  is  then  scaled  so  that 


a  b 


PDFUNS  (a  ,b) =1 . 


In  this  way  a  pseudo-density  function  is  formed.  Finally 
a  uniformly  distributed  random  variable,  r,  [0,1]  is 
generated  and  a  point  (c,d)  is  found  such  that 


c-1 


d-1 


EE,CF“Sa,b*  E 

a=l  b  b=l 

d  d 

PDFUNS 


PDFUNS  ,  <  r 

c ,  b 


(7.8) 


EE 

a=l  b 


,b+  EPDroNSc,b  > r  • 


b=l 


The  Y  associated  with  the  kernel  shape  at  (c,d)  is 
then  used  as  the  next  pixel  in  the  generated  image.  The 
process  is  continued  until  a  full  texture  image  is 
generated  with  the  kernel  window  moving  one  pixel  at  each 
step,  row  by  row. 


In  an  indirect  way,  this  is  equivalent  to  generating 
a  random  variable  having  any  distribution  using  the 
desired  cumulative  distribution  combined  with  a  uniformly 
distributed  random  variable  (which  is  easy  to  generate). 
In  other  words,  uniformly-distributed  deviates  are 
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transformed  to  deviates  having  the  desired  distribution 
using  the  inverse  cumulative  density  function  [45,48,58]. 
This  is  frequently  done  in  simulations. 

7.4  Results 

For  a  kernel  containing  55  pixels  (see  Figure  7.4) 

passed  over  a  128x128  parent  texture  approximately 

7.2xl06  operations  (additions  or  subtractions)  are  needed 

to  get  the  WMSD  image  defined  by  Eq .  (7.4).  Another 

2.6xl05  are  required  to  find  the  next  pixel  according  to 

Eq.  (7.8).  therefore,  to  generate  a  512x512  texture 

12 

requires  1.96x10  (2  trillion)  operations. 

Results  from  texture  synthesis  done  by  this  method 
are  shown  in  Figure  7.5  through  7.15.  The  original 
textures  are  shown  in  Figs.  5.1(a)  through  5.11(a).  Each 
of  these  images  is  512x512  pixels.  A  128x128  section  of 
each  original  (parent)  texture  was  used  for  the 
simulation.  Bark  exhibits  very  large  macro  structure  and 
this  is  lost  in  the  simulation.  A  similar  thing  happens 
with  raffia  as  the  kernel  size  is  smaller  than  the  cell 
size  of  the  original  texture  but  is  not  as  pronounced. 
The  top  part  of  the  bubbles  texture  was  generated  using  a 
128x128  portion  different  than  that  of  the  bottom  part. 
For  this  reason  the  top  20-30%  of  the  texture  looks 
different  from  the  rest.  The  large  number  of  operations 
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Figure  7.3  Perfectly  Correlated  Kernels 
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Figure  7.4  Best-fit  Model  Kernel 


Figure  7.5  Best-fit  Grass  Figure  7.6  Best-fit  Bark 


Figure  7.7  Best-fit  Straw  Figure  7.8  Best-fit  Cloth 
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Figure  7.9  Best-fit  Wool  Figure  7.10  Best-fit 

Leather 


Figure  7.11  Best-fit  Sand  Figure  7.12  Best-fit 

Water 
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Fiqure  7.13  Best-fit  Wood  Figure  7.14  Best-fit 

Raffia 


Figure  7.15  Best-fit  Bubbles 
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makes  this  process  very  time  consuming  even  when 


a 


pipelined  processor  is  dedicated  to  the  task.  About  5.5 
days  of  dedicated  time  on  an  AP120B  were  required  to 
generate  each  texture. 

Although  this  method  is  of  little  practical  use  due 
to  the  computational  complexity  of  the  algorithm  a  few 
points  should  be  made.  With  constantly  increasing 
computer  processing  speeds,  a  simplified  version  of  this 
texture  simulation  method  may  be  implemented  in  the  near 
future.  It  is  even  possible  that  such  computations  could 
be  performed  by  an  array  of  microprocessors.  In  any  case 
such  brute-force  approaches  are  simple  and  could  be  made 
cost-effective  in  the  future. 

The  results  also  indicate  visually  the  amount  of 
texture  information  present  in  a  55  pixel  window  (see 
Figure  7.4)  because  at  each  pixel  generation  step,  the 
next  pixel  in  the  Markov  process  depends  on  only  the 
pixels  in  this  neighborhood. 

Finally,  this  approach  is  admittedly  ad  hoc . 
Numerous  distance  measures  could  replace  the  one  chosen  in 
this  work  and  each  would  give  different  results  that  might 
appear  better  or  worse.  It  is  always  important  that  the 
process  be  random  and  not  merely  copy  the  texture  sample. 
If  the  simulation  region  is  much  larger  than  the  parent 


176 


sample,  a  deterministic  process  will  quickly  generate 
patterns  that  can  easily  be  seen  to  repeat.  In  other 
words,  the  histogram  represented  by  p  ^Vn+1^V1'  *  *  *  ' VN^ 
should  rarely  be  a  delta  function.  A  reduction  in  the 
number  of  computations  could  be  made  if  the  kernel  was 
non-contiguous .  Also,  better  results  could  probably  also 
be  obtained  if  the  kernel  window  was  larger.  The  shape, 
contiguity  and  size  of  the  kernel  in  this  study  was  chosen 
primarily  for  computer  programming  considerations. 

7.5  Conclusions 

The  results  from  this  best-fit  texture  synthesis 
method  are  very  pleasing  but  the  number  of  computations 
required  is  large.  Other  similar  algorithms  could  be 
developed  which  are  simpler  and  could  possibly  produce 
even  better  results.  With  the  decrease  in  computation 
costs  and  the  increase  in  processor  speeds  of  future 
computers,  such  texture  synthesis  methods  could  be 
implemented  in  the  future  without  great  cost. 
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CHAPTER  8 


NON -HOMOGENEOUS  TEXTURES 

8.1  .Introduction 

In  this  chapter,  methods  for  removing  and  introducing 
non-homogeneities  in  texture  images  are  presented. 
Non-homogeneities  in  neighborhood  mean  and  stcindard 
deviation  are  often  removed  previous  to  statistics 
collection  over  a  parent  texture  image  to  improve  the 
accuracy  of  parameter  estimation.  Similar 

non-homogeneities  may  be  added  during  the  texture 
synthesis  process  by  merely  reversing  the  process.  In 
this  way,  synthesized  textures  which  are  homogeneous  may 
be  processed  to  be  non-homogeneous . 

8.2  Removing  Non-Homogeneities 

Prior  to  simulation  attempts,  the  textures  in  this 
study  have  been  preprocessed  by  statistical  differencing 
[42].  This  preprocessing  step  is  described  by 

I(n.,n0)  =  [F(n  ,n  )-F(n.  ,n_)  ]  - — -  + 

2  Ua^.np.oJ  (81) 

[amd+(l-ra)F(n1,n2) ] 
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where  and  represent  desired  mean  and  standard 
deviation.  F  is  the  input  pixel  at  location  (n^,n2>»  row 
n^  and  column  n^  in  the  discrete  digital  image  matrix,  and 
I  is  the  output  pixel  in  the  statistical  differenced 
image.  F(n^,n2)  and  a(n^,n2)  represent  the  mean  and 
standard  deviation  of  the  input  image  at  the  point 
( ni , n 2 ) .  The  variable  A  is  a  gain  factor  that  prevents 
overly  large  output  values  when  a  is  small,  and  a  is  a 
proportionality  constant  controlling  the  extent  to  which 
the  mean  of  the  output  image  is  homogeneous. 

In  our  studies,  mean  and  standard  deviation  factors 
were  computed  in  non-overlapping  16x16  pixel  blocks  values 
are  used  to  compute  the  mean  and  standard  deviation  at 
each  point.  In  our  work  a  =  0.8,  m^  =  128,  85  and 
A  =  6.  These  values  tend  to  induce  local  homogeneity  in 
mean  and  standard  deviation  over  an  image.  Large  amounts 
of  variation,  however,  will  only  be  reduced  and  not 
eliminated  unless  A  is  very  large  and  a  =  1.0. 

A 

Assuming  that  o  does  not  approach  zero,  then  another 
form  of  statistical  differencing  can  be  used.  This  may  be 
written  in  equation  form  as 
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(8.2) 


r  6°d 

I(n1,n2)  =  [F(n1,n2)-F(n1,n2) ]  — - + 

L  6  (n^ ,n2) 

+  [amd+(l-a)F(n1,n2) ] 

Here,  both  a  and  6  are  proportionality  constants  ranging 
from  0.0  to  1.0.  Setting  a.  =  1.0  and  6  =  1.0  returns  an 
output  image  with  precise  desired  mean  and  standard 
deviation.  Setting  a  =  0.0  and  $  =  1.0  causes  the 

standard  deviation  but  not  the  mean  of  the  input  image  to 
change.  Setting  a  =  1.0  and  6  =0.0  causes  the  mean  but 
not  the  standard  deviation  of  the  input  image  to  be 
modified.  Setting  a  =  0.0  and  6  =  0.0  produces  no  change. 
Examples  of  statistical  differencing  are  shown  in 
Figs.  8.1  through  8.4.  The  cork  texture  of  Fig.  8.1  is 
non-stationa ry  in  mean  due  to  shading  differences, 
primarily  at  the  right  edge.  Figure  8.2  shows  the  image 
resulting  from  processing  Fig.  8.1  using  the  statistical 
differencing  algorithm.  The  non-homogeneity  of  mean  is 
removed  and  the  contrast  is  slightly  increased.  An 
original  brick  texture  image  shown  in  Fig.  8.3  has  very 
low  contrast.  After  statistical  differencing,  local 
contrast  is  much  improved  and  the  texture  is  more  apparent 
(see  Fig.  8.4).  Thus  the  statistical  differencing 
algorithm  is  quite  useful  in  eliminating 

non-homogenei t i es . 


(1 


-6) 
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Figure  8.1  Before  Statistical  Figure  8.2  After  Statisti- 
Differencing  cal  Differ¬ 

encing 


Figure  8.3  Before  Statistical  Figure  8.4  After  Statisti- 
Differencing  cal  Differ¬ 

encing 
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8.3  Introducing  Non-Homogeneities 


The  inverse  operation  of  statistical  differencing  can 
be  called  local  moment  modification.  Solving  for  F  in 
terms  of  I  using  Eq.  (8.1)  we  find  the  formula  for  local 
moment  modification  as 


Ac?  (n.  ,n  )+o  F(n1,n,))Aa 

F(ni,n2)  =  - ± — - -  I(nlfn2)+  — ^ — - - 


Aa 


-  (amd+(l-a)F(n1,n2) ] 


Ao  (nlfn2)+ad 


(8.3) 


Using  Eq .  (8.3)  we  can  introduce  non-hornogenei  ti  es  into  a 

simulated  texture  by  generating  an  image  F(n^,n2)  and 

-A 

o(nlfn2)  and  defining  A,  a,  md ,  and  . 

Again,  if  we  assume  that  a( ni>n2^  does  not  approach 
zero,  then  Eq.  (8.2)  can  be  inverted  to  form  another  local 
modification  formula 


F (n^ ,n2)  =  F(nlfn2)  + 


a(n1,n2)  [I (n^ , n2)  -amd~  (1-a) F (n^ , n2 )  ] 


(8.4) 


a (n1,n2) (l-<5)+6c?d 

In  the  process  of  local  moment  modification,  it  is  best  to 
set  m^  and  to  be  equal  to  the  mean  and  standard 
deviation  of  the  homogeneous  texture  I(n^,n2).  Then,  the 
mean  and  standard  deviation  of  the  output  image  (F(n^,n2) 


182 


1 
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1 

1 

! 

will  be  defined  by  the  images  F(n]L,n2)  and  a(nj_,n2). 

These  images  may  be  generated  randomly. 

An  image,  before  and  after  local  moment  modification, 
is  shown  in  Fig.  8.5.  Here  F(n^,n2)  was  assumed  to  be 

A 

ramp-like  and  a(n^,n2)  was  constant.  Many  other  complex 

_ .  /v 

and  random  F(n^,n2)  and  afn^nj)  images  could  be  used  to 
create  different  effects  and  simulate  phenomenon  such  as 
non-homogeneous  lighting  effects. 
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(b)  Synthesis  After  Local 
Moment  Modification 


Cloth 


CHAPTER  9 


TEXTURE  IDENTIFICATION  AND  SEGMENTATION 

9. 1  Introduction 

In  this  chapter  we  examine  various  approaches  to 
texture  segmentation  and  identification  using  the  linear 
model  developed  earlier.  These  methods  employ  covariance 
measures  of  a  texture  over  windows  of  the  region  being 
segmented  or  identified.  These  same  covariance  measures 
were  used  to  estimate  the  linear  model  parameters  which 
were  examined  in  Chapter  3  and  Chapter  5.  In  Chapter  5, 
the  information  content  of  these  measures  was  shown  by 
performing  simulations  of  various  textures  and  therefore, 
based  on  these  results,  we  might  conclude  that  these 
second-order  statistics  would  be  useful  in  the 
segmentation  and  identification  of  textures.  Pictorial 
segmentation  results  are  given  in  this  chapter. 

It  is  generally  agreed  that  a  great  portion  of 
texture  information  is  contained  in  the  second-order 
statistics  of  a  texture.  There  are  notable  exceptions  to 
this  rule  as  was  shown  in  Chapter  2,  however  for  most 
natural  textures,  second-order  statistics  have  proved  to 
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be  sufficient  to  adequately  discriminate  textures  in  most 
applications  [49,51].  In  a  practical  sense,  we  usually 
are  also  prevented  from  gathering  and  analyzing 
higher-order  statistics  because  of  computational 
limitations.  In  fact,  we  can  easily  be  overwhelmed  by 
masses  of  data  arising  from  second-order  statistics. 

The  most  general  second-order  statistics  are  the 

complete  second-order  joint  densities,  or  2-grams,  which 

were  discussed  in  Chapter  2  for  the  one-dimensional 

texture  synthesis  case.  These  measures  may  be  estimated 

by  joint  gray  scale  histograms  taken  over  a  parent 

texture.  For  a  g  grey-level  image,  each  histogram 

2 

requires  g  storage  locations.  But,  as  was  pointed  out  in 

Chapter  2,  there  is  one  such  histogram  for  each  vector 

distance  or  spacing  between  a  pair  of  pixels.  To  compute 

these  histograms  over  all  spacings,  2,...,n  ,  in  two 

2  2 

dimensions  would  result  in  (nr-l)  g  entries.  Even  for 
reasonable  g  and  n  ,  this  could  easily  create  a  data 
expansion  containing  unnecessary  information  rather  than  a 
data  reduction  yielding  measurements  with  discrimination 
value. 

For  these  reasons,  texture  image  data  is  often 

quantized  (to  reduce  g)  and  second-order  measurements 
resulting  in  joint  grey-scale  histograms  are  made  over  a 
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s  in  u  1 1  number  of  pixel  spacinys  (to  reduce  nr)  .  To 
decrease  the  size  of  the  feature  space  further,  various 
functions  of  joint  y rey -scale  his  toy  rams  are  calculated 
and  these  values  are  primarily  used  to  identify  a  texture. 
For  a  sinyle  nistoyram,  as  many  as  25  to  30  functions  have 
been  proposed  [7].  In  spite  of  the  larye  dimensionality 
of  the  feature  space  and  the  problems  with  quantizing  low 
contrast  textures,  this  family  of  texture  features  is  used 
frequently  to  successfully  classify  textures  [7,52]. 

Many  other  identification  and  classification  schemes 
exist  [49]  as  the  discrimination  of  textures  represents 
the  most  important  application  of  texture  analysis. 

In  this  chapter,  we  reduce  the  information  contained 
in  the  joint  yrey-scale  histogram  to  one  single  number, 
the  correlation  coefficient  for  that  particular  pixel 
spacing.  It  is  expected  that  this  large  reduction  will 


cause  a  decrease  in  discrimination  power 

as  the 

size 

and 

information 

content 

of  the  feature 

space 

has 

been 

significantly 

reduced . 

The  purpose  of 

this  exercise 

is 

not  to  develop  a  new,  more  powerful  texture  identifier  but 
merely  to  access  the  information  content  of  the 
correlation  coefficient  values  when  applied  to  the  problem 
of  texture  discrimination.  It  is  already  apparent  from 
the  simulation  results  presented  in  earlier  chapters  that 
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a  great  deal  of  texture  information  may  be  obtained  by 
proper  use  of  these  correlation  coefficients. 

Correlation  measurements  have  been  applied  in  various 
ways  to  the  problem  of  texture  identification  previous  to 
this  study  [2,11,12,14,49,53].  It  has  generally  been 
concluded  that  they  are  of  lesser  value  than  Haralick's 
family  of  functions  on  values  of  the  joint  grey-scale 
histogram  [51]  when  applied  to  the  discrimination  problem. 
In  the  remaining  sections  of  this  chapter  we  will  develop 
two  new  discrimination  methods  utilizing  correlation 
values.  One  is  based  on  the  statistical  test  for  equality 
of  covariance  matrices.  The  other  utilizes  multiple 
statistical  tests  for  the  equality  of  individual 
correlation  coefficients.  Both  show  good  discrimination 
power  but  neither  exceeds  the  quality  of  Haralick's 
measures . 
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pixels  in  length,  was  classified  from  these  measurements 
as  in  Fig .  9.1. 

Second-order  statistics  must  be  measured  over  a 
variety  of  vector  distances  (pixel-pair  spacings)  to  be 
useful  in  texture  discrimination.  In  our  case,  these 
second-order  statistics  are  correlation  values.  There  is 
only  one  correlation  value  for  a  particular  spacing. 

There  are  many  approaches  to  estimating  a  correlation 
value  over  a  window  for  a  particular  pixel-pair  spacing. 
One  involves  passing  a  kernel  of  pixels  over  the  window 
and  taking  a  sample  at  all  points  where  the  kernel  is 
completely  contained  within  the  window  boundaries.  Such 
measurements  would  result  in  a  covariance  matrix  for  the 


kernel  ov 

er  that  window. 

This 

matrix 

can 

be 

used 

to 

identi f y 

a  texture,  as 

will 

be  sho 

wn 

later 

in 

this 

section . 

Another  approach 

to 

estimati 

ng 

a  coi 

'rela 

tion 

value  for  a  particular  pixel-pair  spacing  would  involve 
measurements  over  all  possible  samples  within  the  window 
of  that  spacing.  This  is  equivalent  to  passing  a  kernel 
containing  two  points  over  the  window  for  each  pixel-pair 
spacing.  The  result  of  this  approach  is  a  correlation 
value  for  each  spacing  which  must  be  examined  by  itself 
and  may  not  be  used  to  form  a  correlation  matrix  as  was 
discussed  in  section  3.7.  A  method  for  identifying 
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texture  using  these  individual  correlation  values  will  be 
examined  in  section  9.3. 


Once  we  have  obtained  a  covariance  matrix  by  passing 
a  kernel  of  points  over  a  window,  that  matrix  may  be 
compared  statistically  with  covariance  matrices  from  known 
textures.  Such  a  statistical  test  has  been  derived  for 
testing  the  equality  of  covariance  matrices.  The  maximum 
likelihood  ratio  approach  used  to  derive  this  test  makes 
the  standard  assumptions  of  mul t i-d imens ional  normality 
and  independence  of  samples.  Both  were  used  in  our 
texture  simulation  work.  Details  concerning  the  test  and 
its  derivation  are  given  in  [47,54].  The  statistical  test 
for  two  covariance  matrices  is  given  by 


"l  -  2-3026  dD  1  X2N,N+l)/2 


(9.1) 


where 

D  =  (M1+M2-2)log10|c|-(M1-l)log10|C1|-(M2-l)log10|C2 


C-^  and  C2  are  the  estimated  covariance  matrices  for  each 
g  roup , 

C  =  [ (M1-l)C1+(M2-l)C2l/(M1+M2-2) 

and 

d  -  1  -  [(M^Ty  +  (M^TI  -  (M14-2>][(2"|2-3N-1I/6IN*11.- 

l4  is  the  size  of  the  covariance  matrices  being  tested 
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which  is  equal  to  the  number  of  pixels  in  the  kernel. 
and  M2  are  the  sample  sizes  used  to  estimate  and  C2 . 

| C |  denotes  the  determinant  of  C.  The  test  statistic,  , 
has  an  approximate  chi-square  distribution  with  N(N+l)/2 
degrees  of  freedom  and  approaches  0  as  approaches  C2 . 

Having  derived  the  test  for  equality  of  two 
covariance  matrices  we  must  determine  the  contents  of 
these  matrices.  As  the  number  of  points  in  the  kernel 
increases,  the  size  of  the  covariance  matrices  increases 
leading  to  more  difficult  and  time-consuming  determinant 
calculation  required  in  Eq .  (9.1).  Also,  as  the  spacing 
of  these  pixels  increases,  the  number  of  samples  in  a 
window  decreases.  Finally,  as  more  and  more  points  are 
included  in  the  kernel,  the  amount  of  redundant  and 
overlapping  information  in  the  covariance  matrix  increases 
due  to  redundant  pixel-pair  spacing.  This  was  discussed 
in  section  3.7. 

To  eliminate  this  redundancy,  we  will  consider  the 
problem  briefly  in  one-dimension.  Perfect,  non-redundant 
pixel  spacing  is  possible  only  for  patterns  with  maximum 
range  of  1,  3  or  6  pixels  given  by  the  corresponding 
patterns  XX,  XX-X  and  XX — X-X  as  shown  in  Table  9.1.  For 
these  three  particular  ranges,  the  patterns  shown  are 
constructed  so  that  no  two  pairs  of  X's  are  separated  by 
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redundanc ies 


contain 


if  all  possible  distances  and 


orientations  less  than  or  equal  to  the  maximum  range  are 


to  be  spanned.  Patterns  obtained  using  this  vector 


product  approach  are  heavily  weighted  in  the  horozontal 
and  vertical  directions. 


A  two-dimensional  kernel  containing  16  pixels  which 

spans  a  two-dimensional  range  of  6  is  shown  in  Fig.  9.2. 

Passing  this  kernel  over  a  window  yields  a  covariance 

matrix  of  dimension  16,  which  is  not  an  unreasonable  size 

for  computation  purposes.  The  number  of  data  samples  over 

2 

a  window  of  size  W  is  (W  -6)  .  Even  for  small 

BIG  BIG 

windows  of  size  16,  a  reasonable  sample  size  of  100  may  be 
obtained . 


Once  the  kernel  and  procedure  are  determined,  we 
proceed  with  the  process  of  segmenting  the  textures  to 
test  the  identification  procedure.  There  are  many 
approaches  to  this  problem  which  are  usually  defined  by 
the  particular  case  of  interest.  We  may  or  may  not  have 
prototype  and  parent  texture  data.  We  can  segment, 
cluster-analyze  or  identify.  Notice  that  the  generality 
of  the  test  leaves  a  wide  variety  of  options  open.  We  can 
compare  matrices  in  one  area  of  an  image  with  those  in 
other  areas  and  our  test  statistic  values  will  be  the 
distance  measures  defining  the  closeness  of  the  textured 
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Figure  9.1  Segmentation  Window 


igure  9.2  Two-dimensional  Spanning 
Kernel 


Table  9.1  MINIMAL  SPANNING  SETS 


RANGE 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 


IPTS 

2  XX 

3  XXX 

3  XX -X 

4  XXX-X 

4  XXX— X 

4  XX— X-X 

5  XXXX - X 

5  XXX— X— X 

5  XXX - X— X 

6  XXXX— X - X 

6  XXXX - X - X 

6  XXXX - X - X 

6  XXX - X - X— X 

7  XXXXX - X - X 

7  XXXXX - X - X 

7  XXXX - X - X - X 

7  XXXX - X - X - X 

8  XXXXXX - X - X 

8  XXXXX - X - X - X 

8  XXXXX - X - X - X 

8  XXXXX - X - X - X 

8  XXXX - X - X - X X 

8  XXX - X - X— X— X-X 

9  XXXXXX - X - X - X 

9  XXXXXX - X - - X - X 

9  XXXXX - X - X - X - X 

9  XXXXX - X - X - X - X 

9  XXX - X - X— X— X— XX 

9  XXX - X - X— X— X--X-X 

10  XXXXXX - X - X - X - X 

10  XXXXXX - X - X - X - X 

10  XXXXXX - X - X - X - X 

10  XXXXX - X - X - X - X - X 

10  XXXX - X - X - X - X — X— X 

10  XXX - X - X — X — X — X — X-X 

10  XX-X— X - X - X - X - X - XX 

11  XXXXXXX - X - X - X - X 

11  XXXXXX - X - X - X - X - X 

11  XXXXXX - X - X - X - X - X 

11  XXXXX - X - X - X - X - X - X 

11  XXXX - X - X - X - X - X - X-X 

11  XXXX - X - X - X - X - X— X— X 

11  XX-X— X - X - X - X - X - X XX 

12  XXXXXXX - X - X - X - X - X 

12  XXXXXXX - X - X - X - X - X 

12  XXXXXX - X - X - X - X - X - X 

12  XXXXXX - X - X - X - X - X - X 

12  XXXX - X - X— X - X - X X XX 

12  XXXX - X - XX - X - X - X - X— X 

12  XXXX - X - X - X - X - X - X — X — X 

13  XXXXXXXX - X - X - X - X - X 

13  XXXXXXX - X - X - X - X - X - X 

13  XXXXXXX - X - X - X - X - X - X 

13  XXXXXXX - X - X - X - X - X - X 

13  XXXXXX - X - X - X - X - X - X - X 
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regioi  We  could  compare  a  matrix  with  a  number  of  known 
prototypes  and  classify  the  unknown  region  according  to 
the  test  statistic  for  each  comparison.  Finally,  we  may 
merely  wish  to  locate  a  particular  texture  in  an  image  and 
thus  we  would  compare  matrices  with  only  one  prototype. 

This  later  approach  was  selected  in  this  chapter  as 
the  visual  results  are  simple  to  display  and  analyze.  We 
will  compare  the  matrix  measured  over  a  large  window  with 
one  obtained  over  the  complete  parent  texture  prototype 
and  assign  the  pixels  in  the  small  window  accordingly  (see 
Fig.  9.1).  Figure  9.3  shows  the  test  composite  image 
used.  We  will  attempt  to  identify  the  texture  sand  in 
that  region.  The  texture  sand  was  chosen  as  the 

simulation  results  of  this  texture  were  in  some  ways  very 
poor  when  examined  critically  on  a  high  resolution  device. 
Therefore,  it  represents  a  worst-case  example. 

Figure  9.4  shows  the  segmentation  results  when 

W  =  W  =  16.  The  pixel  brightness  are  linearly 

mapped  test  statistic  values,  which  are  supposed  to  be 

chi-square-distributed.  Improvement  is  made  when 

W  =32  and  W  =16  as  is  seen  in  Fig. 9.5. 

BIG  SMALL 

Figure  9.6  shows  much  improved  results  when  W  =  48  and 

*  BIG 

w„  =  16.  The  difference  is  due  to  the  availability  of 

SMALL  1 

more  texture  information  in  the  larger  window  in  the  form 


Figure  9.5  Segmentation  Using  Figure  9.6  Segmentation  Using 
Covariance  Matrix,  Covariance  Matrix, 


W 


BIG 


32 


WBIG  =  48 
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of  increased  sample  size  (100  to  1764). 

As  a  final  note,  it  should  be  understood  that  this 
chosen  texture  kernel  (see  Figure  9.2)  totally  ignores 
information  obtained  earlier  concerning  the 
interdependence  of  pixels  in  the  generating  kernel  as 
expressed  in  the  linear  model.  It  is  very  possible  that 
better  segmentation  results  could  be  obtained  by  using  the 
kernel  shape  which  best  fits  each  texture  as  this 
particular  pattern  of  points  was  found  to  be  most 
significant.  The  linear  model  used  for  simulation  of  each 
texture  might  even  be  used  itself  as  the  d is tr ibut ions  of 
the  parameter  estimates  are  known  if  certain  assumptions 
are  made  [24,39].  But  neither  of  these  approaches  have 
been  tried  even  though  improvements  are  expected. 

9.3  Segmentation  Using  Individual  Correlation 
Coefficients 

Individual  correlation  estimates  for  particular 
pixel-pair  spacings  may  be  made  by  taking  measurements 
over  all  possible  samples  within  the  window  containing 
that  spacing.  This  approach  utilizes  more  information 
than  the  fixed  kernel  method  as  it  includes  measurements 
near  all  edges  of  the  window  (thus,  the  sample  size  is 
increased) .  The  result  is  a  single  coefficient  for  each 
pixel-pair  spacing.  If  the  spacing  between  each  pair  is 
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unique  then  the  set  of  coefficients  will  be  non- redundant . 
The  covariance  matrices  used  to  segment  textures  in  the 
last  section  may  contain  redundant  information  if  the 
pixel-pair  spacings  in  the  kernel  are  repeated. 

A  statistical  test  for  comparison  of  two  correlation 
coefficients  has  been  derived  for  the  bivariate  normal 
distribution  and  is 


where 

U2  =  (Z1-Z) 2 (M1-3)+ (Z2-Z) 2 (M2-3) 

and 

Z ^  =  arctanh  (r^) 

Z2  =  arctanh  (r2) 

Z3  =  [(M1-3)Z1+(M2-3)Z2l/(M1+M2-6) 
r^  and  r 2  are  the  computed  correlation  coefficients  from 

the  two  groups.  This  test  is  derived  using  the  fact  that 

is  approximately  distributed  N(Z-^:  arctanh  ( P-^)  , 

-1 

(M1-3)  )  where  is  the  true  correlation  coefficient  of 

the  population  (see  section  3.5).  For  additional  details 
see  [39].  Again,  to  derive  the  tests  assumptions  of 
normality  and  independent  samples  are  made. 

As  in  the  previous  section,  the  test  statistics  may 
be  used  to  segment  or  identify  textures.  However, 
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difficulty  arises  as  our  test  must  be  performed  for  each 
measured  correlation  value.  Thus  a  number  of  U2  values 
are  obtained  if  multiple  pixel-pair  spacings  are  used.  At 
this  point  these  values  can  be  combined  in  numerous  ways. 
As  each  value  is  chi-square  distributed  we  can  compute 

the  probability,  prT  ,  that  a  random  variable,  which  has 

u2 

that  chi-square  distribution,  is  greater  than  or  equal  to 
U2  as  shown  in  Fig.  9.7.  For  U2  =  0  (r-^  =  r  2}  this 
probability  is  1.0.  A  product  of  the  probabilities 
associated  with  each  U2  was  used  to  indicate  the  overall 
fit  of  one  set  of  correlation  coefficients  to  another  set. 

Finally,  it  should  be  noted  that  tests  for  equality 
of  correlation  coefficients  will  not  detect  differences  in 
mean  and  variance  over  the  window.  The  correlation 
coefficient  specifically  removes  these  effects.  As  a 
result,  it  may  be  advantageous  to  include  tests  for 
equality  of  means  and  variances  into  the  comparison 
process.  For  equality  of  means  the  test  statistic  is 


are  the  calculated  variances  for  the  two  groups  being 
compared.  This  test  statistic  is  a  value  of  a  random 
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variable  having  the  t-distr ibut ion  with  M^+M2~2  degrees  of 
freedom.  For  the  equality  of  variances  the  test  statistic 
is 


t-,  = 


1 

Sf  1  Mj-l.Hj-i 


(9.4) 


This  test  statistic  is  a  value  of  random  variable  having 

an  F-d  istr  ibut  ion  with  M^-l  and  f^-l  degrees  of  freedom. 

The  derivation  of  both  of  these  tests  requires  the 

assumption  of  normality  for  the  two  groups  and  an 

independent  sample  set.  The  known  distributions  of  the 

test  statistics  t  and  f^  was  used  to  obtain  the 

probabilities,  p  and  p_  that  a  random  variable  having 

1 

that  distribution  would  be  greater  than  or  equal  to  the 
calculated  values  of  1 1 1 1  and  MAXff-^l/f^)  or  less  than 
—  1 1  ^ |  and  MIN  ( f  ^  ,  1/f  ^)  (see  Figs.  9.8  and  9.9). 


Having  obtained  the  set  of  probabilities,  p  , 

U2 


and 


p  and  pf  corresponding  to  the  set  of  tests  for  the 
tl  1 

equality  of  correlation  coefficients  and  the  tests  for 
equality  of  mean  and  variance  of  the  window  a  single  test 
statistic 


'y  ;log  (MAX  (pT1  ,10  6))+log  Pt  +log  p 


+  log  p 

1  rl 

(9.8) 

fit  of  these 

tests. 

to  eliminate 

seal  inj 
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problems  which  would  occur  when  taking  the  product  of  many 


small  numbers.  The  p7  's  are  bounded  below  to  reduce  the 

2i 

effect  of  single,  noisy  correlation  coefficients.  These 

output  test  statistics  obtained  by  comparing  a  windowed 

region  of  Fig.  9.3  with  statistics  gathered  from  the 

parent  prototype  texture  sand  were  linearly  scaled  to 

produce  the  intensity  images  shown  in  Figs.  9.11  to  9.13. 

The  results  for  W  =W  =16  are  shown  in  Fig.  9.11. 

BIG  SMALL 

Improved  results  for  WBIG  =  32  and  WSMALL  =  are  shown 

in  Fig.  9.12.  Results  for  W'  =  48  and  W  =16  are 

J  BIG  SMALL 

shown  in  Fig .  9.13. 


9.4  Conclusions 


In  implementing  the  two  procedures  detailed  in  this 
chapter,  large  values  for  the  test  statistics  U and  U2 
were  obtained  even  when  comparing  the  matrices  or 
correlation  values  measured  from  an  extracted  portion  of  a 
parent  texture  with  those  obtained  from  the  entire  parent 
texture.  There  are  two  basic  reasons  for  this.  First, 
the  assumptions  of  normality  are  probably  incorrect. 
There  is  little  reason  to  believe  that  the  distribution  of 
pixels  in  an  image  is  truly  multi-variate  normal. 
Secondly,  the  samples  used  to  estimate  the  covariance 
matrix  and  the  correlation  coefficients  are  not  random  and 
independent.  They  are  strongly  related  by  their  spatial 
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Figure  9.10  Input  Composite 
Image 


Figure  9.11  Segmentation 
Using  Indivi¬ 
dual  Correla¬ 
tion  Values, 


Figure  9.12 


Segmentation 
Using  Indivi¬ 
dual  Correla¬ 
tion  Values, 


32 


Figure  9.13 


Segmentation 
Using  Indivi¬ 
dual  Correla¬ 
tion  Values, 


48 


positions  in  the  image  and  are  often  adjacent.  The  sample 
size  may  be  adjusted  to  reflect  this  fact  but  no  work  in 
this  area  was  done. 

In  most  cases,  statistical  tests  may  be  considered  to 
be  robust  -  not  weakened  significantly  -  when  assumptions 
used  to  derive  them  are  violated.  This  is  probably  not 
true  of  the  test  for  equality  of  two  covariance  matrices 
which  has  been  called  "frail  at  best"  [47],  For  these 
reasons,  the  methods  presented  may  be  considered 
innovative  but  not  necessarily  statistically  sound. 

In  spite  of  this  fact,  the  results  show  that  the 
method  using  the  test  for  equality  of  covariance  matrices 
was  superior  to  the  method  involving  multiple  tests  of 
individual  correlation  coefficients.  This  could  be  due  to 
the  method  used  to  combine  the  multiple  tests. 

The  pictorial  results  of  this  chapter  indicate  the 
usefulness  of  correlation  values  in  texture 
identification.  The  methods  are  not  intended  to  be 
improvements  over  existing  segmentation  identification 
techniques.  An  adequate  number  of  discrimination 
techniques  have  been  proposed  by  researchers  already.  The 
discussion  is  intended  to  apply  a  specific  texture 
synthesis  model  which  was  based  on  second-order  statistics 

to  the  texture  identification  problem  and  this  was  done 
successfully. 
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CHAPTER  10 


SUMMARY 

10.1  Introduction 

In  this  chapter,  the  results  of  each  of  the  texture 
synthesis  methods  presented  in  this  thesis  are  discussed 
and  compared.  This  summary  will  clarify  the  methods  and 
analyze  the  favorable  and  unfavorable  characteristics  of 
each . 

10.2  Tabulation  and  Discussion  of  All  Models 

A  complete  synopsis  of  the  synthesis  methods 
presented  in  this  thesis  is  effectively  contained  in 
Table  10.1.  Listed  in  this  table  are  the  eleven  methods 
of  texture  generation  and  simulation  presented  in 
Chapter  2  through  Chapter  7.  The  first  three  columns  of 
the  table  state  the  location  of  the  text  and  figures 
(output)  associated  with  each  and  also  the  key  figure  or 
equation  number  which  identifies  the  method  in  a  simple, 
straightforward  manner. 

The  next  three  columns  are  used  to  evaluate  the 
complexity  of  statistics  collection,  statistics  storage 


16.M*  20  i.y.  10  brut-fore,  woro-h 


and  the  generating  process  as  presented  in  Fig.  1.1.  The 
computational  requirements  are  approximate  and  are 
indicated  in  CPU  time  of  a  DEC  KL-10.  Naturally,  these 
numbers  are  relative  to  the  processor  used.  The  storage 
requirements  are  given  in  terms  of  full-word  processor 
locations  needed  for  statistics  storage.  In  some  cases, 
this  storage  could  be  reduced  by  packing  more  than  one 
number  (especially  integers)  into  one  32-  or  36-bit  word 
but  that  was  not  done  here.  For  example,  four  8-bit 
integer  values  will  fit  into  one  32-  or  36-bit  word. 

The  seventh  column  provides  a  relative  measure  of  the 
quality  of  the  texture  simulation  on  a  zero  to  ten  scale. 
A  value  of  5  indicates  a  good  or  reasonable  simulation. 
As  synthesis  evaluation  is  a  nebulous  process  so  is  the 
assignment  of  relative  merit.  The  assessment  of  results 
is  internal  to  this  thesis  as  there  is  little  synthesis 
work  in  the  general  literature. 

The  last  column  of  the  table  contains  important 
general  comments  on  each  method. 

The  one-dimensional  binary  generation  method  of 
Chapter  2  permitted  study  of  visual  response  to  changes  in 
the  probability  distribution  of  texture.  Although  the 
method  was  not  useful  for  natural  texture  simulation,  it 
did  lead  to  the  models  of  later  chapters. 
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The  N-gram  model  of  Chapter  3  was  an  extension  of  the 
Chapter  2  model  applied  to  two-dimensional  natural 
textures.  The  results  were  very  good  even  with  the  severe 
constraint  imposed  by  the  upper  limit  on  the  number  of 
pixels  allowed  in  the  generation  kernel.  With  an  increase 
in  the  complexity  of  the  collection  process  this  model  was 
extended  (see  section  3.9)  to  allow  both  a  larger  number 
of  pixels  to  be  in  the  kernel  and  an  increase  in  kernel 
range.  The  extension  produced  better  simulations  of 
structured  textures.  Finally,  in  section  3.10,  the  binary 
linear  model,  which  was  used  to  determine  the  contents 
(shape)  of  the  generation  kernel,  was  used  to  generate 
binary  textures.  The  textures  generated  using  this  model 
were  nearly  equal  in  quality  to  those  of  the  more  complex 
and  storage-consuming  N-gram  model.  The  N-gram  model  of 
Chapter  3  uses  a  generation  kernel  whose  contents  (shape) 
depends  on  the  linear  model.  Therefore,  the  number  of 
computations  required  in  the  statistics  collection  portion 
of  the  N-gram  model  necessarily  includes  computations  of 
the  linear  model.  However,  in  some  cases,  points  which 
lie  far  from  the  kernel  eye  can  be  neglected  in  the  N-gram 
model  as  only  the  best  few  are  used  due  to  storage 
limitations.  On  the  other  hand,  such  points  should  be 
included  in  the  linear  model  therefore  a  larger 
neighborhood  surrounding  the  kernel  eye  should  be  used  in 
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the  estimation  of  the  linear  model. 

Realizing  the  power  of  second-order  statistics,  a 
method  of  reproducing  Nth  order  statistics  using  algebraic 
reconstruction  was  presented  in  Chapter  4.  The  approach 
proved  to  be  academic  as  the  number  of  iterations  leading 
to  a  solution  yielding  adequate  synthesis  results  was  very 
large  and  required  much  storage  and  computation.  The 
simulations  were  also  slightly  less  appealing  than  the 
methods  of  Chapter  3. 

In  Chapter  5,  the  linear  autoregressive  model  of 
Chapter  3  was  applied  to  256-gray-level  imagery.  The 
results  were  good  considering  the  vast  reduction  of 
information  caused  by  the  statistics  collection  process. 
Slightly  better  results  were  obtained  by  allowing  the 
model  to  contain  cross  terms  but  the  resulting  complexity 
suggests  that  the  change  in  texture  quality  is  not  worth 
the  added  effort  and  computational  expense.  Using 
non-gauss  ian ,  non-stationary  noise  in  the  model  (see 
section  5.7)  produced  markedly  better  results  but  with  a 
requirement  of  slightly  increased  storage. 

The  ski p-generate  method  of  Chapter  5  may  be  used  to 
improve  the  simulation  of  textures  having  a  coarse 
structure.  The  model  produces  results  equal  in  quality  to 
the  linear  autoregressive  model  of  Chapter  5  while 
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requiring  fewer  computations  in  the  collection  process. 
The  piecewise  linear  autoregressive  model  presented  later 
in  Chapter  6  promises  an  analytically  superior  fit  but 
produced  results  which  were  not  appealing  enough  to 
warrant  additional  effort. 

The  field-definition  model  presented  at  the  end  of 
Chapter  6  is  useful  when  generating  textures  composed  of 
subtextures.  The  idea  of  defining  fields  for  texture 
generation  is  a  powerful  approach  for  both  statistics 
collection  and  texture  synthesis  but  the  slow  response  of 
the  autoregressive  model  to  boundaries  produced  results 
less  appealing  than  expected. 

The  best-fit  model  of  Chapter  7  represents  a 
brute-force  approach  to  texture  synthesis.  Though 
computationally  demanding,  the  final  results  show  that 
excellent  texture  simulations  can  be  generated  using 
complete  statistics  from  a  relatively  small  neighborhood. 
The  problems  with  a  small  neighborhood  are  seen  in  the 
simulation  of  regular  textures  such  as  raffia  where  the 
size  of  the  primitives  in  the  texture  is  much  greater  than 
the  window  used  in  the  best-fit  calculation. 

Combining  the  method  of  choosing  a  kernel  shape  using 
a  linear  model  (discussed  in  Chapter  3)  with  the 
ski p-generate  and  best-fit  models  would  result  in  a  very 
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powerful,  but  extremely  computationally-demanding  texture 
synthesis  method.  The  effort  required  to  generate 
textures  by  such  a  complex  combination  of  methods  would 
not  be  required  in  many  cases  where  simpler  models  will 
produce  adequate  results.  The  complexity  of  the  model 
depends  on  the  texture  being  simulated. 

The  ideas  presented  in  Chapter  8  allow  the 
introduction  and  removal  of  texture  non-homogeneities.  In 
Chapter  9,  texture  identification  methods  using  the 
statistics  employed  in  model  parameter  calculation  of 
earlier  chapters  were  proposed.  Although  they  do  not 
discriminate  textures  as  well  as  the  methods  of  other 
researchers  they  do  illustrate  the  application  of 
synthesis  models  to  texture  identification  problems. 

10.3  Suggestions  for  Future  Study 

Many  approaches  to  texture  analysis  have  been  carried 
out  in  the  frequency  domain  by  using  transform  and 
filtering  techniques.  Texture  synthesis  could  be  carried 
out  in  such  a  domain  or  on  frequency-filtered  image.  For 
example,  numerous  textures  could  be  generated  in 
non-overlapping  frequency  planes  and  then  added  together 
to  obtain  a  final  texture  synthesis.  However,  each  of 
these  planes  is  probably  interdependent  and  a  simple 
generation  with  summation  is  probably  not  possible. 
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Another  future  approach  might  employ  to  a  greater 
extent  those  statistics  useful  for  texture  identification 
and  discrimination.  In  most  cases,  these  measurements  are 
not  readily  suitable  to  a  synthesis  process  but  with 
careful  study,  many  could  possibly  be  used  in  such  a 
manner.  Still,  there  is  little  evidence  thus  far  to 
indicate  that  statistics  useful  for  texture  identification 
will  be  useful  for  texture  synthesis. 

An  area  which  deserves  immediate  attention  involves 
preprocessing  of  texture  by  convolution.  Noise  filtering 
and  smoothing  could  be  useful  in  improving  model  parameter 
estimates.  Also,  a  deconvolution  processing  of  images 
synthesized  to  resemble  convolved  textures  might  produce 
excellent  simulation  results. 

Along  this  line,  another  synthesis  method  should  be 
considered.  With  any  model  there  is  always  unexplained 
variance  which  the  model  fails  to  account  for  when  it  is 
applied  to  the  original  parent  texture.  Ideally,  this 
unexplained  variance  would  be  merely  noise  but  this  is 
rarely  the  case  in  practice.  It  is  possible  to  develop 
additional  models  which  could  be  used  to  explain  (or 
simulate)  this  prev iously-unexpla ined  variance.  Combining 
these  new  models  with  the  original  model  would  result  in 
an  improved  synthesis  method. 
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ti  .ally,  texture  is  often  produced  by  the  shading 
effect  of  a  light  source  on  a  three-dimensional  object. 


Therefore,  to  generate  a  texture,  a  three-dimensional 
relief  could  be  formed  then  shaded  using  current, 
commonly-used  graphics  techniques.  Such  an  approach  could 
be  worthwhile  in  some  cases  but  in  most,  synthesizing 
three-dimensional  relief  (which  is  the  same  as  generating 
height  information  over  a  two-dimensional  grid)  is 
equivalent  to  generating  intensity  information  over  a 
two-dimensional  grid. 

10.4  Conclusion 

Many  natural  textures  are  generated  using  a  variety 
of  methods  presented  in  this  thesis.  The  quality  of  the 
natural  texture  simulations  depends  on  the  amount  of 
computation  and  storage  used  in  each  process.  Many 
textures  were  adequately  simulated  using  simple  models 
thus  providing  a  potentially  great  data  compression  for 
many  applications.  Others  required  more  extensive 
computation  to  synthesize  visually  pleasing  results. 
Thus,  as  might  be  expected,  the  success  of  any  synthesis 
method  is  highly  dependent  on  the  texture  itself.  When 
examining  the  results  of  any  method  the  characteristics  of 
both  the  model  and  the  textures  used  must  be  considered'. 
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It  would  be  unwise  to  believe  that  all  textures  could 
be  generated  using  any  single  approach,  especially  one 
which  promises  to  compress  texture  information  to  a 
handful  of  numbers.  Yet  this  is  precisely  what  has  been 
attempted  in  the  texture  synthesis  work  of  th.s  thesis. 

It  is  important  to  note  the  power  and  complexity  of 
each  synthesis  method  of  this  thesis.  Many  textures  can 
be  simulated  well  using  simple  models  such  as  the 
autoregressive  model  if  the  model  is  carefully 
constructed.  Improvements  in  texture  simulation  were  made 
by  modifying  these  models  and  allowing  them  to  become  more 
complex  and  use  more  information  in  the  generation 
process.  Other  textures  require  more  complex  models  such 
as  the  best-fit  model  of  Chapter  7.  The  shortcomings  of 
each  method  will  constantly  indicate  where  future  work  can 
be  done. 
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APPENDIX  A 


GENERATING  BINARY  TEXTURE  PAIRS  POSSESSING 
EQUAL  SECOND-ORDER  STATISTICS 

Texture  pairs  may  be  generated  which  have  equal 
second-order  statistics  but  are  visually  d  iscr iminable. 
Let  G  (V,,  V  2>  v3'  V4^  represent  the  probability  of 
generating  a  0  after  the  pixels  ^j'V2'V3'V4  a^on<3  a  l*ne 
in  a  one-dimensional  texture  (a)  and  G i' V2' V3' V4 ^ 
textures  ( b) .  (This  is  a  slight  change  of  notation  from 
Chapter  2  as  (a)  and  (b)  will  be  texture  number  indices  in 
this  Appendix.)  Define 

V-  =  | l“Vi |  ( A . 1 ) 

where  V^e{  0,1). 

The  restrictions  used  to  generate  textures  with  equal 
second-order  statistics,  pa(vg»vj)  =  pb^Vl,Vj^'  *n  Chapter 
2  may  be  stated  as 

Ga(Vl'V2'V3'V  =  Ga(Vl'V2'V3'V4)  (A. 2) 

WWV  =  1~Ga<Vl'V2-V3'V  <fl-3) 
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and 


vww  = 


Gb  *V1'V2'V3,VV 


(A. 4) 


Equations  (A. 2)  and  (A. 3)  must  also  hold  for  texture 
(b) .  Combining  Eq .  (2.14)  and  (2.24)  for  a  4-gram  system, 
it  can  be  shown  that 


Pa(V 

Pa(V 


1 


1 


'V3'V 


'V3'V 


P 

a 


(V1'V2'V3,V4) 

(v1,v2,v',v4)  . 


(A. 5) 

(A. 6) 


Combining  the  above  restrictions  yields  two  further  rules, 

V  V' 

[V5+(-l)  5Ga(V1,V2,V3,V4)]=[V5+(-l)  5Gb(V1,V2,V3,V4) ](A-7) 

V  v' 

[V5+(-l)  5Ga(V1,V2,V3,V4)]=[V'+(-l)  5Ga(V1,V2,V3,V4) ] 

v' 

=  [V'+(-l)  5Ga(V;[,V2,V3,V4)](A*8) 
These  restrictions  yield  textures  having  the  equalities 


pa(°)=  Pb(0)  =  Pa(l)  =  Pb(l)  •  (A. 9) 

As  no  closed-form  solutions  to  Eq.  (2.14)  and  Eq.  (2.24) 
are  easily  obtained  the  proofs  verifying  Eq.  (A. 5)  and 
(A. 9)  are  very  complex.  However,  these  properties  may  be 
illustrated  by  generating  numerous  examples  where 
Eq.  (A. 2),  Eq.  (A. 3)  and  (A. 4)  hold. 
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Assuming  the  above  conditions,  a  proof  then  follows 
to  show  that 
restrictions  hold: 


Pa(VL,Vj)  =  P^fV^fVj)  if  the  above 


Pa(VVj>  -  E  Pa<VV2 . Vl'V 

V. 

1 


i^i  /  j 


“IX'VWvJ 

*XIPb(Vl-V2'V3'V4)J  'VW)'5cblV*'Vj'V2’Vl11 


n  tVk+(-l)  "Ga(V  V  V  V  )] 
k=5 

j 
n 

k=5 


=y^Ph(vI,v^,v3,v4)  . 


(A. 10) 


M 


/ 

VV-_o+(-l) 


v 


3k+2 


kI5V"3b-2 


Gb(V3k-2,V3k-l,V3k,V3k-l) 1 


[V3k-3+(-1)  3Gb(V3k-l'V3k,V3k-l,V3k-2) ] 


|V3k+4+  <-X  >  <V3k '  V3k-1  'Vik-2 '  V3k+3 

v: 


»> 


^  [vj+l-H+<-1)  ^+l"^b(vj-t-3'Vj-e-2'Vj-«-l'Vj-l) ' 
where  r  =  M0D(j-l,3).  It  is  very  important  to  note  that 
the  variables  in  the  product  expression  match 
successively.  That  is,  V^k-l  anc^  V3k+2  d  ^  have  the  prime 
notation  in  each  sub-expression.  At  this  point,  there  are 
three  possible  cases: 


0 


Case  1 ,  r  = 


In  this  case,  the  proof  is  completed  by  a  series  of 

change  of  variables  as  summing  over  \A  is  the  same  as 

summing  over  V' .  The  remainder  of  the  proof  becomes 

V, 


Ed 

Pb(Vl'V2'V3,V4) ‘  n  tvk+(-1)  Gb(Vk-4,Vk-3'Vk-2,Vk-l 
k=l  ,  k 


)] 

(A. 11) 


-  w 


Case  2 ,  r  =  1 


In  this  case,  the  remainder  of  the  proof  becomes 

j-1  vk 

Pb(Vl'V2'v3'V4}  n  ^vk+(-1)  Gb(Vk-4'Vk-3'Vk-2'Vk-l) 1  * 


k=l 


V  . 


•IV1’  3Gb(Vj-4'Vj-3’Vj-2-Vj-l 


)] 


j-1 


(A. 12) 


Pb(Vl'V2'V3'V4)lrn1IVk+(“1)  Gb(Vk-4'Vk-3'Vk-2'Vk-l)] 

k=l 

•IVj4W)\fVj.4.V].3.V).J,V  l(l  . 


Thus,  in  this  case  Pa(vl'vj)  =  pb*vl'vj)'  This  implies 
that  Pa(01r0j)  =  Pb(01flj);  Pa^O^lj)  =  pb^°i'°j^;  Pa^1r°j^ 
=  Pb(lrlj)»  Pa(lrlj)  =  Phav0.).Vie  know  Pa(Orlj)  -  Pa(Vj) 
and  PfcdL^j)  and  P  (0 )  =  P  (1 )  =  0 . 5.  So  Pa(VlfVj) 

=  Pb(Vi,Vj)  =  0.25  for  all  Vx  and  Vj  . 
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Case  3 ,  r_  =  2_ 

In  this  case,  as  in  the  case  where  r  =  0, 
remainder  of  the  proof  is  straightforward,  requiring  only 
a  change  of  variables. 
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APPENDIX  B 

GOODNESS-OF-FIT  TEST 
FOR  N-GRAMS 

We  may  think  of  the  one-dimensional  binary  texture 
generation  process  as  an  experiment  which  should  yield 
certain  N-grams  given  certain  generation  parameters. 
These  N-grams  or  Nth-order  densities  may  be  determined 
analytically  by  combining  Eq .  (2.14)  and  Eq.  (2.24)  of 
Chapter  2.  However,  as  with  any  Monte-Carlo  simulation, 
the  analytic  parameters  may  not  agree  exactly  with  the 
statistics  or  estimates  of  those  parameters  based  on 
observations  from  the  simulation  or  experiment. 
Statistical  tests  may  be  used  to  determine  whether  the 
statistics  match  the  parameters  to  the  extent  expected 
given  a  given  random  sample. 

The  most  common  statistical  test  for  this  purpose  is 
the  "goodness-of-f it"  test  involving  the  chi-square 
distribution.  This  test  is  used  when  we  want  to  compare 
an  observed  distribution  with  the  corresponding  values  of 
a  theoretical  distribution.  The  test  is  based  on  the 
statistic 
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2N  ..  .2 


(B.l) 


The  and  the  e^  are  the  observed  and  expected 

frequencies  of  a  certain  N-gram  pattern  in  our  experiment. 

The  sampling  distribution  of  this  statistic  is 

N-l 

approximately  the  chi-square  distribution  with  2 
degrees  of  freedom.  (The  constraints  on  the  system  of 
N-grams  yield  this  number  of  degrees  of  freedom.)  A 
chi-square  distribution  with  n  degrees  of  freedom  is  given 
by 


2  (V2) ^2  (n-2)  _ y 2/2 

f(X  )  =  -  l  X  X  ( B .  2 ) 

2n/2r (J) 

The  approximation  is  close  provided  e  >5. 

As  a  first  step  in  the  testing  process,  the  complete 
set  of  N-grams,  P  ^Vl  '  V2  '  •  •  •  , VN  )  i-5  computed  using 
Eq .  (2.14)  and  Eq .  (2.24).  Their  corresponding 

statistics,  P (V^ , V2 , . . . , VN)  are  computed  from  the 
generated  texture  by  counting  the  number  of  occurrences  of 
each  pattern  (V^,...,VN)  and  dividing  by  the  total  sample 

A 

size,  M.  Then,  fi  =  P (V^, . . . , VN) *M  and 

e^  =  P  (V^,  . .  .  ,V^)  *M.  For  large  N,  the  number  of  degrees 
of  freedom  is  greater  than  30  and  a  normal  distribution 
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is 


qui  te 


accurate . 


The 


expression 


approximation 


- 


/2N-1  is  approximately  normally  distributed  as  the 


standard  normal  distribution.  The  tables  listed  in  [18] 
may  also  be  used  for  large  N.  Additional  details 


concerning  the  goodness-of-f it  test  may  be  found  in 
[45,63-65] . 


For  each  of  the  textures  in  Chapter  2,  this  test  was 
applied  to  each  texture  half  and  the  hypothesis  that  the 
expected  and  measured  N-grams  are  statistically  equal  was 
accepted  at  a  a =  0.01  confidence  level  for  N  =  1,...,10. 
Actually,  these  tests  are  not  independent  in  a  statistical 
sense  as  the  information  content  for  each  N  overlaps.  But 
this  poses  no  violation  of  assumptions  of  the  statistical 
test.  However,  all  possible  patterns  in  the  generated 
sample  may  be  a  violation  of  the  independence  of  samples 
assumption.  Still,  a  low  chi-square  value  indicates  a 
good  fit  of  the  data  to  predicted  parameter. 


224 


APPENDIX  C 


PSEUDO-RANDOM  VARIATE  GENERATION 

Two  algorithms  were  used  in  this  thesis  to  generate 
uniformly-distributed  pseudo- random  numbers.  Both  are 
modifications  of  the  multiplicative  congruence  method 
[48,66].  These  uniformly-distributed  variates  may  be 
transformed  to  normally-distributed  deviates  using  the 
inverse  normal  probability  distribution  function. 

The  first  algorithm  is  the  Lehner  Pseudo-Random 
Number  Generator  [67].  The  general  form  of  this  generator 
is  - 

Xn+1  =  KVM0D  231_1)  (C.l) 

where 

K  =  14 2 9 (MOD  231-1)  =  63036001610 

The  resulting  number  is  basically  a  31-bit  pattern  which 
is  in  the  form  of  a  floating-point  number  in  the  range 
(0,1).  The  algorithm  can  be  written  to  avoid  division  by 
231-1  [68  1. 

The  second  algorithm  used  to  generate  uniformly 
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distributed  variates  is 

Xn+1  =  ,5  Xn(M0D  231-1)  '  (C.2, 

_  O  I 

The  resulting  integer  is  multiplied  by  2  to  convert  it 
to  a  floating-point  number.  This  generator  is  reported  in 
Refs.  [69-71]. 

These  uniform  deviates  may  be  transformed  to 
normally-distributed  deviates  if  desired.  this  is 
accomplished  by  computing  the  inverse  of  the  integral 


2/2 

dt 

( C .  3 ) 

✓H  J 

~oo 

The 

algor 

'ithm  for 

accompl ish ing 

this 

is 

descr 

ibed 

in  Ref. 

[57] 

and 

is  briefly  outlined  here. 

The 

basic  interval  (0,1)  is 

d  iv 

ided 

into 

4 

segments . 

In 

each 

segment 

the  inverse 

of 

the 

Gauss 

ian 

integ  ral , 

invgauss(X),  or  an  integral  of  a  similar  form  is 
approximated  by  a  minimax  rational  function  [55,72].  The 
approximations  for  the  final  3  segments  (comprising  the 
interval  Xe{(0, 0.075)  (0.925,1)})  are  functions  of  the 
transformed  variable 

W  =  V-log  e (1-X2)  .  (C . 4 ) 
This  transformation  of  the  variable  improves  the 
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efficiency  and  stability  of  the  approximation  [SO],  The 
rational  functions  in  these  intervals  are  of  the  form 

(C^W+C-W  +C3W 

cn+  - T - 3 

dQ+d1W+d2W  +W 

For  the  remainder  of  the  (0,1)  interval, 
xe{(0,075,  0.925)  },  the  function  is  of  the  form 

invgauss(X)  = 

( C  .  6  ) 

(Z+Z- (b0+a1.Z2/(b1+Z2+a2/(b2+Z2+a3/(b3+Z2) ) ) ) )  *SGN(Z) 

where  z  =  1 1  —2  x  | .  The  constants  a^,  bj_,  c  ^  and  d^  have 

specific  values  found  by. solving  the  minimax  problem  and 
are  given  in  Ref.  [57]. 


)■ 


(C .  5 ) 
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